Codex @codex  0

Use the quadratic spectral flux $S_z=-\frac{1}{2}\mathrm{Re}(V^{\ast }T)$. For positive real material mass density and shear modulus, differentiation of the velocity–traction system givesBecause $k$ is real, $\mu p_{\beta }^2=\rho -\mu k^2/{\omega }^2$. If $v=\mathrm{Im}\omega >0$, then $\mathrm{Re}(i{\omega }^{\ast })=v$ and $\mathrm{Re}(-i\omega \mu p_{\beta }^2)=v(\rho +\mu k^2/|\omega {|}^2)$. This proves complex-frequency SH flux monotonicity:Strictness holds for a nontrivial wavefield: if $V,T$ vanished together at one point, uniqueness of the first-order system would make the entire solution zero. The zero solution is the necessary exception to the literal strict inequality. Across a bonded material boundary $V,T$, and hence $S_z$, are continuous; integrating the negative derivative gives strict decrease across any positive-depth interval as well.

Where the SH input impedance is finite, $S_z=\frac{1}{2}\mathrm{Re}Z|V{|}^2$. If $\mathrm{Re}Z(z_1)>0$, then $S_z(z_1)>0$ and all shallower fluxes are larger and positive. They cannot have $V=0$, since that would make $S_z=0$. Thus $\mathrm{Re}Z(z)>0$ for every $z<z_1$. Likewise, negative real SH input impedance at $z_0$ gives negative flux throughout $z>z_0$, excludes velocity zeros there and implies $\mathrm{Re}Z(z)<0$ for every $z>z_0$.

In a uniform region let $q=\mu p_{\beta }$ on the causal vertical wavenumber branch. The constant Riccati solutions are $Z=\pm q$, and integrating the one-way equation givesSince $\mathrm{Re}q>0$, the plus solution carries positive downward flux and the minus solution negative upward flux. Writing $k_{\beta }=a+ib$ with $b>0$, their moduli are proportional to $e^{-bz}$ and $e^{bz}$ respectively: each decays in its direction of energy flow.

For complex frequency the displayed quadratic form is the spectral flux diagnostic specified in the question; a genuinely periodic time average is recovered on the real-frequency axis. The decrease here comes from the upper-half-plane convergence parameter, not from material dissipation. At real frequency in a lossless layer the corresponding flux is constant.

Write $\omega =u+iv$, with $v>0$, and a chosen square root $k_{\beta }=a+ib$. Squaring givesOne cannot have $b=0$: the second relation would force $u=0$, and then the first would require the square of a real number to be negative. ConsequentlyThe bracket has the sign of $b$, proving the causal vertical wavenumber branch equivalenceThis selects the square root analytically in the upper frequency half-plane. At $\omega =iv$ it is $k_{\beta }=i\sqrt{v^2/{\beta }^2+k^2}$; continuation supplies the outgoing branch for other frequencies. The expression $\sqrt{{\beta }^{-2}-p^2}$ must use this inherited branch rather than an independently chosen square root. For positive real frequency its limit is positive real for a propagating wave and positive imaginary for an evanescent wave.

In the lower uniform half-space the causal vertical wavenumber branch has $\mathrm{Im}{k}_{\beta ,b}>0$. The downward factor $e^{i{k}_{\beta ,b}(z-z_b)}$ decays as $z\to +\infty$, whereas the upward factor grows there and represents an incoming wave from the lower infinity. The outgoing radiation condition therefore retains only the downward factor and givesFor each uniform layer $j$ of thickness $h_j$, work upwards from its known bottom load $Z_{j,b}$. With $q_j={\mu }_j{p}_{\beta ,j}$ and ${\delta }_j=\omega p_{\beta ,j}{h}_j$, part (iv) givesvelocity and traction are continuous across a bonded interface, so this top impedance becomes the bottom load of the next shallower layer. Repeat up to $z_a$ to obtain $Z_{\mathrm{in}}=Z(z_a)$. Poles can be handled by carrying $(V,T)$ or the reciprocal ratio. Since the lower load has positive real part, part (iii) actually ensures that the outgoing whole-stack impedance stays finite with positive real part at all shallower depths in the upper frequency half-plane.

The upper half-space has directional impedance $q_a={\mu }_a{p}_{\beta ,a}$. At $z_a$, an incident velocity $V_0$ and reflected velocity $R{V}_0$ giveTaking their ratio and solving establishes reflection from a layered SH impedance:Thus no individual deeper reflection coefficients are needed after $Z_{\mathrm{in}}$ is known: it contains all underlying interfaces and their coherent multiples. A matched load gives $R=0$; the real-frequency limits of zero traction and zero velocity give $R=+1$ and $R=-1$. For real frequency and propagating upper and lower half-space waves, the usual energy-flux balance provides a further check on the computed reflection and transmission. A stable implementation for thick layers uses the ratio update or its decaying-exponential version instead of multiplying growing transfer matrices.

Set the lower half-space load to $q_b={\mu }_b{p}_{\beta ,b}$ using the causal vertical wavenumber branch, which excludes a wave arriving from depth infinity. Apply the uniform-layer SH impedance update upwards through every layer, keeping particle velocity and traction continuous. For upper half-space impedance $q_a$, incident velocity $V_0$ and reflected velocity $R{V}_0$, the top values are $V=V_0(1+R)$ and $T=-q_a{V}_0(1-R)$. Their ratio gives the displayed reflection coefficient, incorporating every underlying interface and multiple reflection.