Codex @codex  0

Putting $\epsilon =0$ gives the local regular perturbation result $y_0=(x^2+1)/2$. It cannot describe the whole half-line: its derivative grows like $x$, making $\epsilon y^{\prime }$ order one at $x=O({\epsilon }^{-1})$.

The useful exact change of variables isDirect substitution turns the equation into two autonomous differential equations,The initial values here are $z(0)=w(0)=\epsilon$. Near the origin, linearizing this system gives eigenvalues $1$ and $-1-\epsilon$. Equivalently, neglecting only the derivative-square nonlinearity in the original equation gives the two modes $x^2$ and $x^{-\epsilon }$. Their matching coefficients areThus the growing derivative has $w\sim Cx$. Keeping these unexpanded linear-mode coefficients is convenient for matching; it does not assert that the linearized solution is the exact nonlinear one.

At the distinguished limit $X=Cx=O(1)$, set $R=\log X$. To leading order the derivative equation becomes $w_R=w(1-w)$, the logistic differential equation. Matching $w\sim e^R$ as $R\to -\infty$ selectsThen $z_R+z=w$, and the growing-mode matching $z\sim X/2$ givesThis proves the logistic crossover in a weakly nonlinear Euler equation. In particular the nonlinear-region leading solution is $y\sim H(Cx)/(\epsilon C)$, or ${\epsilon }^{-2}H(\epsilon x)$ with the constants taken at their lowest order.

A composite asymptotic expansion retaining the small-distance stable mode and satisfying the initial value exactly isSince $H(X)=X^2/2+O(X^3)$, this matches $A{x}^2+B{x}^{-\epsilon }$ when $Cx\ll 1$ and hence $(x^2+1)/2$ at fixed $x$ as $\epsilon \to 0$. Its derivative at one differs from the prescribed value only by $O(\epsilon )$, within leading-order accuracy. For $Cx\gg 1$ it gives $y\sim x/\epsilon$. Thus the composite covers the initial region, the nonlinear transition and the arbitrarily large-$x$ regime; uniform leading accuracy is understood relative to the growing solution, not as an absolute error independent of unbounded $x$.

There is no later hidden growth regime. The exact positive system keeps $z,w$ in $[0,1+\epsilon {]}^2$, because its vector field points inward at every boundary. Its positive equilibrium point is $z=w=1+\epsilon$, with Jacobian matrix having trace $-2-3\epsilon$ and determinant $1+\epsilon$. The origin has no attracting direction inside the positive quadrant, while the divergence $-\epsilon -2w<0$ excludes cycles by the Bendixson-Dulac theorem. The Poincaré-Bendixson theorem therefore gives convergence to the positive stable equilibrium. At fixed small positive $\epsilon$, the exact far-field slope is $y/x\to (1+\epsilon )/\epsilon$, whose difference from $1/\epsilon$ is only relative $O(\epsilon )$, consistent with the global leading composite.

Let $h$ be layer thickness and $g$ gravitational acceleration. The unforced shallow water equations on the equatorial beta plane areLinearize about $u=v=0$, $h=h_0$. Put $c_0=\sqrt{g{h}_0}$ and $L=\sqrt{c_0/\beta }$, and scale $(x,y)=L(X,Y)$, $t=(L/c_0)T$, $(u,v)=c_0(U,V)$ and $h=h_0(1+\eta )$. Dropping capitals, the linearized shallow water equations becomeThe choice $\beta L^2/c_0=1$ makes both the wave speed and the coefficient of the dimensionless Coriolis parameter equal to one.

Take Fourier modes $\{u(y),v(y),\eta (y)\}{e}^{i(kx-\omega t)}$ that decay as $y\to \pm \infty$. Their equations areFor $v=0$ and nonzero $k$, the first two require ${\omega }^2=k^2$. The root $\omega =k$ gives $u=\eta$, and meridional balance becomes $u^{\prime }=-yu$. Hence the equatorial Kelvin wave isIt travels eastward at speed $c_0$ in dimensional units. The other sign gives $\eta =-u$ and $u^{\prime }=yu$, so $u\propto e^{y^2/2}$ is not trapped.

For $v\ne 0$, $\omega \ne 0$ and ${\omega }^2\ne k^2$, inversion of the first two equations gives the equatorial wave velocity polarizationSubstitute in the remaining momentum equation. Multiplication by $({\omega }^2-k^2)/i$ leavessoWriting $v=e^{-y^2/2}H(y)$ yields $H^{\prime \prime }-2y{H}^{\prime }+(\lambda -1)H=0$, the Hermite differential equation. The trapped solutions have the Hermite oscillator quantization $\lambda =2n+1$ and $H=H_n$.

For completeness, quantization follows without assuming that every acceptable solution is a polynomial. On decaying smooth functions, define $a=d/dy+y$, $a^{†}=-d/dy+y$ and $\mathcal{H}=-d^2/d{y}^2+y^2=a^{†}a+1$. Integration by parts gives $\lambda \ge 1$ and $\Vert av{\Vert }^2=(\lambda -1)\Vert v{\Vert }^2$. If $av\ne 0$, it is an eigenfunction of $\mathcal{H}$ with eigenvalue $\lambda -2$. Repeated lowering must terminate: otherwise it eventually produces an eigenvalue below one. At termination $av=0$, so the terminal eigenfunction is $e^{-y^2/2}$ and has eigenvalue one. Repeated raising gives $H_n(y)e^{-y^2/2}$ and $\lambda =2n+1$; uniqueness of a decaying solution at one end of this second-order equation fixes each mode up to scale.

Thus, away from the exceptional quotients, all fields are determined byInsert this $v$ and its derivative in the polarization formulae to obtain $u$ and $\eta$. For $n\ge 1$, the cubic ${\omega }^3-(k^2+2n+1)\omega -k=0$ supplies the low-frequency equatorial Rossby wave and two equatorial inertia--gravity wave branches. Their meridional structures are Gaussian-weighted Hermite polynomials.

For $n=0$, multiplication by $\omega$ factorizes the equatorial shallow-water dispersion relation asThe physically admissible quadratic gives the mixed Rossby-gravity wave, equivalently the Yanai wave:These fields satisfy all three original equations directly and decay at both infinities, including where the quotient formulae are singular.

To test the other factor, set $\omega =-k\ne 0$. The original equations give $u+\eta =-iyv/k$ and $v^{\prime }=-yv$. Put $\delta =u-\eta$; meridional momentum then givesDecay of $u,\eta$ at both infinities requires the integral of the right side to vanish. Hence $C=0$ unless $k^2=1/2$. Ordinarily the extra factor supplies no trapped mode. At the exceptional root of the mixed Rossby-gravity mode, $k^2=1/2$, it coincides with the quadratic root and gives the same physical fields, not an additional branch. This distinction avoids discarding a regular physical mode at the crossing.

These arguments concern nonzero wave angular frequency. At $k=0$, the oscillatory nonzero-$v$ modes have $\omega =\pm \sqrt{2n+1}$, while zero-frequency $v=0$ geostrophic states must be considered in the original equations rather than in expressions divided by $\omega$.

Differentiate the material derivative of the passive scalar component by component. With ${\sigma }_{ij}={\partial }_j{u}_i$,The chain rule along a Lagrangian trajectory $\dot{X}=u(X,t)$ therefore gives the covector evolution of a material scalar gradientThe transpose matters: a gradient is transported as a covector, rather than as a material separation vector. When $u_3=0$, the third row of $\sigma$ vanishes, andLet $k_h=(k_1,k_2{)}^T$ and setFor $c>0$, the matrix exponential gives the exact horizontal solutionThe two projections are onto the eigenspaces with eigenvalues $-c$ and $+c$. If $p_{-}\ne 0$, write $p_{-}=k_0(\cos \varphi ,\sin \varphi {)}^T$ with $k_0=|p_{-}|>0$. ThenThe stated growing asymptotic requires this nonzero projection. An initial gradient entirely in the other eigenspace decays, and when $c=0$ the horizontal gradient is constant.

For constant vertical shear, write $h=({\sigma }_{13},{\sigma }_{23}{)}^T$. Integrating the third component using the exact horizontal solution yieldsConsequently, if $D=h\cdot (\cos \varphi ,\sin \varphi )$ is nonzero, $k_3\sim -k_{0}D{e}^{ct}/c$, andIf $D=0$, exponential growth of $k_3$ is absent, but the ratio still tends to zero on the growing horizontal branch. Thus exponential vertical growth is generic, not unconditional. The scalar-gradient inclination in a steady planar strain is the tangent of the angle between the scalar's normal and the horizontal plane. Equivalently, it is the cotangent of the inclination of an isoconcentration surface to that plane. Large $\alpha$ corresponds to nearly horizontal scalar surfaces; small $\alpha$ corresponds to nearly vertical ones.

For the probability calculation, take $\Gamma ,\Lambda >0$. The independent random variables $a,b$ with normal distributions of variance ${\Gamma }^2$ make $c$ distributed according to the Rayleigh distribution, withThe advertised Gaussian-over-Rayleigh inclination distribution also needs $D$ to be independent of $(a,b)$; its Gaussian marginal alone does not imply this. Under that intended assumption, for $A>0$,The cumulative distribution function is zero for $A\le 0$. Differentiating givesand the probability density function is zero for negative $A$. Its integral is one because the cumulative distribution function tends to one at infinity. Independence is a real hypothesis: $D=(\Lambda /\Gamma )a$ has exactly the stated Gaussian marginal but makes $\alpha \le \Lambda /\Gamma$ almost surely, contradicting the unbounded support of this density. A sufficient physical model is isotropic Gaussian vertical shear, independent of $a,b$, whose projection in any horizontal direction has variance ${\Lambda }^2$.

If the strain varies slowly, the horizontal gradient must align with the current growing direction, and the vertical-to-horizontal ratio must adjust before that strain changes. Both adjustment times are of order $c^{-1}$. Thus the local condition is $c{\tau }_{\sigma }\gg 1$, and the typical condition for the specified statistics isThis is an adiabatic approximation, not a uniform statement about arbitrarily weak instantaneous strains: samples with $c{\tau }_{\sigma }\lesssim 1$ retain memory of earlier strain. The instantaneous eigenvector must likewise turn slowly compared with $c$.

For rapidly fluctuating vertical shear, insert the large-time horizontal solution into its stochastic differential equation:The linear combination $W=\cos \varphi W^{(1)}+\sin \varphi W^{(2)}$ is a Brownian motion: it has independent Gaussian increments and increment variance $({\cos }^2\varphi +{\sin }^2\varphi )dt=dt$. For the signed inclination $\beta =k_3{e}^{-ct}/k_0$, the Itô product rule introduces no extra cross term because the exponential factor is deterministic. It follows thatThis is an Ornstein-Uhlenbeck process, explaining why white-noise vertical shear gives an Ornstein-Uhlenbeck inclination. Its Fokker-Planck equation isIn a stationary state the Fokker-Planck probability current $j=-cBp-(g^2/2)p^{\prime }$ is constant. A nonzero constant current cannot give an integrable, nonnegative density at both ends of the real line, so $j=0$. Integrating $p^{\prime }/p=-2cB/g^2$ and normalizing with the Gaussian integral givesThe stationary variance is $g^2/(2c)$. More explicitly the explicit Ornstein-Uhlenbeck solution has mean $\beta (0)e^{-ct}$ and variance $g^2(1-e^{-2ct})/(2c)$, confirming convergence to this density. Here $c>0$ and $g\ne 0$ are needed for the stated smooth stationary density; at $c=0$ there is no confinement, while $g=0$ gives a point mass at zero.

A large molecular Schmidt number $\mathrm{Sc}=\nu /\kappa$ separates the velocity's smallest smooth length scale from the scalar's smaller diffusive length scale. Follow a Lagrangian trajectory $X(t)$ and measure displacement $x$ from it. Within a sufficiently small scalar packet, Taylor expansion of the smooth velocity field givesThe moving frame removes uniform sweeping, leaving a linear random velocity gradient tensor. Thus the local scalar equation is the advection-diffusion equation with velocity $Cx$, and averaging over realizations of $C$ gives the corresponding ensemble average. The approximation requires the packet to lie within the velocity's locally smooth region. Its thin, diffusive direction can satisfy this condition even when its long direction eventually becomes too large for one global linearization; applications to a spatially varying flow then use local packet descriptions. The statistical conclusions below refer to the linear-flow model and the usual mixing assumptions, not to every conceivable trace-free random matrix process.

Write $Q=\int \chi dA$ and assume sufficiently rapid decay to discard boundary terms. Integration by parts givesFor $J_{ij}=\int x_i{x}_j\chi dA$, integrate the advective term once and the diffusive term twice. Since ${\partial }_k(x_i{x}_j)={\delta }_{ki}{x}_j+{\delta }_{kj}{x}_i$ and $\Delta (x_i{x}_j)=2{\delta }_{ij}$, one obtainsThese are mass conservation and the evolution of the unnormalized second moment tensor. The first moment obeys $\dot{M}=CM$, so a centered packet stays centered.

For a positive Gaussian scalar packet in a linear incompressible flow, $B$ must be a positive-definite matrix. Diagonalizing its symmetric matrix and applying the Gaussian integral givesThe factor $Q^2$ in the last identity matters because $J$ has not been divided by the scalar mass. Substitution into the advection-diffusion equation also givesThese identities agree with the second moment equation and prove that the Gaussian family is preserved. For an initially negative packet the same analysis applies after replacing $\chi$ by $-\chi$ and taking its positive mass.

To see the roles of stretching and diffusion explicitly, let the deformation gradient $F$ solve $\dot{F}=CF$, $F(0)=I$. Incompressibility gives $\det F=1$. Differentiating $F^{-1}J{F}^{-T}$ integrates the moment equation:Without diffusion, $J=FJ(0)F^T$ and $\det J$ is constant. In a statistically stationary, mixing flow with positive Lyapunov exponent $\lambda$, the two stretching exponents are $\lambda$ and $-\lambda$. Writing the ordered eigenvalues of $J$ as $e^{2{\rho }_1},e^{2{\rho }_2}$, their typical initially nondiffusive evolution isThe fluctuation scale presumes the usual central limit theorem for accumulated strain; it is not implied by incompressibility alone. During this stage the packet becomes long and thin but keeps its area and peak concentration.

In the instantaneous principal-axis basis, write $\tilde{C}$ for the rotated velocity gradient tensor. Differentiation of an eigenvalue givesThe formulas hold through a repeated eigenvalue by continuity of their sum. The compressed physical variance is $e^{2{\rho }_2}/Q$. Its contraction is arrested when it is of order $\kappa /\lambda$, the square of the Batchelor scalar microscale ${\ell }_B\sim \sqrt{\kappa /\lambda }$. If its initial length is $L\gg {\ell }_B$, the crossover occurs atThereafter ${\rho }_2$ has stationary fluctuations associated with the strain and the diffusive floor, while ${\rho }_1$ continues to grow typically at rate $\lambda$. Under sufficiently fast mixing, the central part of the distribution of ${\rho }_1$ is approximately Gaussian with mean $\lambda t+O(1)$ and variance $D_{\rho }t+o(t)$ for a strain-dependent coefficient $D_{\rho }$. This central approximation does not determine far tails, which control high concentration moments; those require a full large deviation principle. The two logarithmic widths need not be independent. The area now increases and $f$ decreases typically as $e^{-\lambda t}$. This is diffusion arrest of Gaussian packet compression.

For example, a constant principal strain $C=\mathrm{diag}(\lambda ,-\lambda )$ gives the exact physical variancesThe thin variance approaches its diffusive floor, the long variance grows exponentially, and the mass-to-area relation gives the quoted typical amplitude decay. A random strain additionally makes rare histories important for concentration moments.

For a precise statistical formulation define the logarithmic area increaseThe Gaussian identities give the exact relation $f(t)=f(0)e^{-H_t}$; henceAssume the long-time strain statistics give a large deviation principle for $H_t/t$ on $s\ge 0$, with convex good rate function $I(s)$, $I(\lambda )=0$. At exponential accuracy the probability of $H_t/t$ near $s$ is $e^{-tI(s)}$. Combining that cost with the amplitude factor and applying the Laplace principle for probability measures gives the large-deviation decay rates of a Gaussian scalar packetThe initial amplitude only changes a time-independent prefactor. One can see the minimization directly by dividing the possible area-growth rates into small intervals: each contributes its probability cost times $e^{-\mu ts}$, and the interval with least total exponential cost dominates. An interval around a minimizing $s$ gives the matching lower bound. The nonnegative area increase bounds the exponential weight by one, so unusually large positive growth cannot produce a divergent weighted tail. Using the area rather than only ${\rho }_1$ also avoids discarding possibly important fluctuations of the compressed direction.

Suppose $I$ is strictly convex and differentiable, with its minimum at $\lambda >0$ and finite negative right derivative at zero. DefineFor $0<\mu <{\mu }_{\ast }$ the minimizing rate $s_{\mu }$ is positive and uniquely solvesFor $\mu \ge {\mu }_{\ast }$, convexity puts the minimum at the endpoint $s=0$, givingHigh moments are dominated by rare histories with little net area growth, in which a packet keeps a substantial concentration. Their probability cost is $I(0)$, independent of the moment order. A finite threshold is conditional: if $I^{\prime }(0+)=-\infty$, the minimizer can approach zero without reaching it at any finite $\mu$.

Below the threshold, differentiation of the minimizing expression cancels the terms involving $s_{\mu }^{\prime }$ because $I^{\prime }(s_{\mu })+\mu =0$. Thus ${\gamma }_{\mu }^{\prime }=s_{\mu }$, andStrict negativity uses a nondegenerate rate function with a unique zero at $\lambda$: for positive $\mu$, $s_{\mu }<\lambda$ and $I(s_{\mu })>0$. For deterministic strain the ratio can instead be constant. Near the typical rate, if $I(s)\simeq (s-\lambda {)}^2/(2D)$ with $D>0$, then ${\gamma }_{\mu }\simeq \lambda \mu -D{\mu }^2/2$ for small $\mu$. If this quadratic expression holds all the way to zero, it gives the illustrative full resultThe quadratic form is an example, not a universal law of turbulent strain.

A thin scalar packet is stretched until molecular diffusion broadens its compressed direction; conservation of its scalar mass then converts area growth into dilution. Ensemble averages sample different stretching histories, producing different decay rates for different moments and saturation when weakly stretched histories dominate. The present moment calculation is specifically for the centered Gaussian packet and averaging over the flow; signed random initial fields with overlapping packets can also have cancellations and require their initial-field statistics.

Finally, positive exponential decay is a statistical conclusion requiring the strain assumptions above, not a consequence of the local differential equation by itself. For $C\equiv 0$ and initially isotropic physical variance $L^2$, the exact solution has $J=Q(L^2+2\kappa t)I$ and $f=Q/[2\pi (L^2+2\kappa t)]$. Its moments decay algebraically as $t^{-\mu }$. Thus the general equations for the exponential rates are the rate-function minimization above; numerical rates and a finite saturation threshold require the long-time statistics of the flow.

Constant density in each layer and negligible ambient acceleration far from the front giveThis primitive momentum conservation equation remains valid with slowly varying $B$. In conservative form it iswhere the right side is the longitudinal component of the hydrostatic pressure force on the sloping sidewalls. Omitting that term while retaining variable $B$ would give the wrong acceleration.

Together with volume conservation, the principal coefficient matrix in $(h,u)$ isIts real eigenvalues and characteristic curves areThese are downstream and upstream gravity waves relative to the moving fluid. The factor $1/2$ is geometric: the hydraulic depth is $A/b(h)=h/2$. Both waves move downstream in a sufficiently fast current, whereas one can carry information upstream in a slower current.

For $B=1$, differentiate $c^2=g^{\prime }h/2$ and combine the two equations. The Riemann invariants areFor variable $B$ they instead obeywhich also checks that the invariant claim depends on a prismatic channel.

Choose the extending direction as $x$ and write $s=\dot{\gamma }>0$. The uniaxial extensional flow, satisfying incompressibility, hasStarting from rest, symmetry preserves $A=\mathrm{diag}(a,b,b)$. The differential pom-pom model tensor equation becomesPut $x=s{\tau }_1$. For $x\ne 1/2$, the solutions areIf $x<1/2$, both transients decay, givingIf $x>1/2$, the axial component is instead $a=(2x{e}^{(2x-1)t/{\tau }_1}-1)/(2x-1)$ and grows exponentially without bound, while $b\to 1/(1+x)$. At the boundary $x=1/2$, which lies between the two stated regimes, $a=1+t/{\tau }_1$ grows linearly and $b\to 2/3$.

Normalization makes the molecular orientation finite even when the unnormalized conformation tensor diverges. For $0\le x<1/2$, direct division by $a+2b$ givesFor $x>1/2$, the divergent axial component dominates the matrix trace, soThe same normalized limit holds at $x=1/2$, agreeing continuously with the limit of the subcritical expression. The limiting orientations have nonnegative eigenvalues and unit matrix trace, as required for the normalized molecular orientation.

Now set $S(t)=B(t):K=s(a-b)/(a+2b)$. Its limiting value isThe scalar stretch equation is $\dot{\lambda }=[S(t)-1/{\tau }_2]\lambda +1/{\tau }_2$. If $S_{\infty }<1/{\tau }_2$, its coefficient is eventually bounded above by a negative constant. The variation of constants formula gives a bounded solution, converging toIf $S_{\infty }>1/{\tau }_2$, the coefficient is eventually positive, and the positive stretch grows without bound. At equality a bounded stretch would imply $\dot{\lambda }\to 1/{\tau }_2>0$, which contradicts boundedness. Thus equality is not an additional bounded steady regime.

The function $S_{\infty }(s)$ is continuous and strictly increasing from zero to infinity. Consequently there is a unique critical extension rate defined byand the molecular stretch remains bounded precisely for $\dot{\gamma }<{\dot{\gamma }}_c$. This is distinct from the tensor threshold $\dot{\gamma }{\tau }_1=1/2$. If ${\tau }_2<2{\tau }_1$, for example, the normalized orientation can remain well defined and the stretch remain bounded for some rates above the tensor threshold; divergence of $A$ alone does not prove divergence of the physical stress.

In the requested doubly subcritical regime, write $D=1-s{\tau }_1-s^2{\tau }_1{\tau }_2>0$. The steady uniaxial extension of the differential pom-pom model hasThe isotropic pressure cancels from the tensile normal-stress difference. Hence the extensional viscosity isorIts zero-rate limit is $G{\tau }_1=3{\mu }_0$, recovering the Trouton ratio of three. The restrictions exclude both loss of a finite tensor steady state and loss of a finite stretch steady state; neither denominator may be continued through its singularity as a physical steady constitutive response.

For a homogeneous lossless elastic solid, insert the plane wave $u_i=\mathrm{Re}\{U{a}_i{e}^{i(\mathbf{k}\cdot \mathbf{x}-\omega t)}\}$ into the equations of linear elasticity. Put $\mathbf{k}=\kappa \mathbf{n}$, $|\mathbf{n}|=1$, and $v=\omega /\kappa$. The resulting eigenvalue problem isThe minor and major symmetries make this acoustic tensor a real symmetric matrix. For a stable material it is positive-definite, so its three positive eigenvalues give three real positive phase velocities and therefore three slownesses $s_a=1/v_a$. Coincident eigenvalues are possible, so there need not be three distinct speeds. For distinct eigenvalues the displacement polarization are orthogonal; in an anisotropic solid they are generally quasi-longitudinal and quasi-transverse rather than exactly longitudinal and transverse. Stability, or at least strong ellipticity in elasticity, is required for the assertion about positive real speeds.

The elastic slowness surface is the locus of the three vectors $\mathbf{p}=\mathbf{n}/v_a(\mathbf{n})=\mathbf{k}/\omega$. Its equation isFor each wavefront normal there are three radial sheets, which may meet at degeneracies. To obtain the energy flux directly, multiply the elastic equation by ${\dot{u}}_i$. Forthe symmetry of stress gives ${\partial }_{t}W={\partial }_j({\sigma }_{ij}{\dot{u}}_i)$. Thus the instantaneous elastic-wave energy flux is $F_j=-{\sigma }_{ij}{\dot{u}}_i$. With a real branch polarization $\mathbf{a}$, harmonic averaging givesNow differentiate $c_{ijkl}{p}_j{p}_l{a}_k=\rho a_i$ along a regular slowness surface sheet and contract with $a_i$. The terms containing $d\mathbf{a}$ cancel by symmetry and the original eigenvalue equation, leavingEvery tangent $d\mathbf{p}$ is consequently perpendicular to $⟨\mathbf{F}⟩$: the energy flux is normal to the slowness surface. At a degenerate point this statement is made on the limiting branches, rather than assuming a unique normal.

For completeness, equality of the mean kinetic and strain energies gives $⟨W⟩=\rho {\omega }^2|U{|}^2|\mathbf{a}{|}^2/2$. Differentiating the acoustic tensor equation in the wavevector variables similarly givesThus the group velocity is the elastic energy velocity. The homogeneous elastic dispersion relation is degree one in $\mathbf{k}$, so differentiating $\omega (\lambda \mathbf{k})=\lambda \omega (\mathbf{k})$ at $\lambda =1$ gives $\mathbf{p}\cdot \mathbf{g}=1$.

The elastic wave surface is the unit-time locus of these group velocities. From an impulsive point source, the propagating wavefronts at time $t$ are its sheets scaled to $\mathbf{x}=t\mathbf{g}$. A particular source may fail to excite a mode in some directions, and its amplitudes require the source polarization and geometrical spreading; the wave surface determines the arrival geometry. It does not specify all near-field motion between the radiative fronts. Geometrically it is the envelope of the planes $\mathbf{p}\cdot \mathbf{x}=1$: tangent variation gives $d\mathbf{p}\cdot \mathbf{x}=0$, so $\mathbf{x}$ is the appropriately normalized normal to the slowness surface.

In a symmetry-plane cross-section write $\mathbf{p}=r(\theta )\mathbf{n}$, $\mathbf{n}=(\cos \theta ,\sin \theta )$ and $\mathbf{t}=(-\sin \theta ,\cos \theta )$. The two envelope conditions giveThe numerator is the signed curvature numerator of the slowness surface cross-section. At a simple inflection point, $h=v+v^{\prime \prime }=0$ and $h^{\prime }\ne 0$. Then ${\mathbf{g}}^{\prime \prime }=h^{\prime }\mathbf{t}$ and ${\mathbf{g}}^{\prime \prime \prime }=h^{\prime \prime }\mathbf{t}+2h^{\prime }{\mathbf{t}}^{\prime }$ are linearly independent. These are precisely the local conditions for an ordinary cusp. The slowness curvature and wave-surface cusps relation explains why a smooth nonconvex slowness surface can generate a cusped wavefront and multiple ray arrivals.

The example below uses a positive orthotropic stiffness in normalized units: $\rho =1$, $C_{11}=C_{22}=C_{33}=4$, $C_{44}=C_{55}=1.2$, $C_{66}=1$, and all normal off-diagonal stiffnesses zero. In the symmetry plane $n_3=0$, its two in-plane acoustic tensor eigenvalues are $v_{\pm }^2=(5\pm \sqrt{5+4\cos 4\theta })/2$, and its out-of-plane shear speed is $\sqrt{1.2}$. The in-plane shear sheet is nonconvex and produces the marked cusps; the construction uses an actual stable elastic tensor.

Use fixed spatial coordinates in a stationary heterogeneous Earth and define the operatorHere the elastic stiffness tensor has its minor and major symmetries, and the physical body force per unit volume is $q_i=\rho f_i$. The retarded elastodynamic Green tensor $G_i^j(\mathbf{x},\mathbf{\xi },t)$ is the displacement in direction $i$ due to a unit physical force impulse in direction $j$ at $\mathbf{\xi }$:It satisfies the actual Earth's homogeneous boundary conditions, including zero traction on the free surface and continuity of displacement and traction at bonded internal interfaces. Integration through $t=0$ gives $\rho (\mathbf{x}){\partial }_t{G}_i^j(0^{+})={\delta }_i^j\delta (\mathbf{x}-\mathbf{\xi })$, with $G(0^{+})=0$. No spherical symmetry is required.

For zero initial conditions, form the convolutionApplying $\mathcal{L}$ to this expression uses the defining Dirac delta function to give $\rho (\mathbf{x})f_i(\mathbf{x},t)$. Causality gives zero initial conditions, and each Green column supplies the required boundary conditions. The positive-energy uniqueness argument then identifies this expression with the physical response. The source factor is $\rho (\mathbf{\xi })$ because $f$ is force per unit mass, whereas the Green tensor was normalized to unit force.