Radiation Forces in Scattering Media as Virtual Work

May 20, 2026

1. Motivation

The Generalized Wigner–Smith (GWS) framework connects incident wavefronts to radiation forces on objects embedded in complex scattering environments. Existing derivations rely on field-theoretic variational identities within either acoustics or electromagnetism separately. Here we show that the central result — the identification of the GWS matrix expectation value with a generalized radiation force — follows directly and in full generality from the principle of virtual work, with the port phase shifts serving as the natural energy bookkeeping variable. This approach (i) unifies acoustic and electromagnetic cases without invoking Maxwell’s stress tensor or acoustic radiation pressure explicitly, (ii) makes the independence of virtual displacements transparent, and (iii) extends naturally to lossy media, where a richer force decomposition emerges.

2. Lossless Case

2.1 Setup

Consider a lossless scattering region coupled to the environment through $$N$$ ports. Incident and outgoing wave amplitudes at port n are $$a_n$$ and $$b_n$$ , normalized so that $$|a_n|^2$$ and $$|b_n|^2$$ equal the corresponding power flows. The scattering matrix $\mathbf{S}$ is unitary ( $\mathbf{S}^\dagger \mathbf{S} = \mathbf{I}$ ) and relates them via $\mathbf{b} = \mathbf{S}\,\mathbf{a}$ .

The GWS matrix associated with a parameter $\alpha$ (e.g., the position of an embedded object) is

\mathbf{Q}_\alpha = i\,\mathbf{S}^\dagger\,\partial_\alpha \mathbf{S}. \tag{1}

2.2 Virtual Work Principle

For a set of objects with positions $\{\mathbf{r}_i\}$ , the principle of virtual work states

\delta W = -\delta U, \tag{2}

where the virtual work done by radiation forces is

\delta W = \sum_i \mathbf{F}_i \cdot \delta\mathbf{r}_i. \tag{3}

Notice the $\delta \mathbf{r}_i$ are completely independent from each other, which means the principle should yield $$N$$ independent equations if there are no additional geometry constraints.

2.3 Energy Bookkeeping via Port Phases

Write each outgoing amplitude in polar form: $b_n = \sqrt{P_n}\,e^{i\phi_n},$ where $$P_n = |b_n|^2$$ is the outgoing power and $\phi_n$ its phase. A virtual displacement $\delta\mathbf{r}_i$ perturbs the scattering medium, changing both powers and phases at the ports. For a single-frequency (harmonic) excitation, the energy stored in the scattering region associated with port n is related to the group delay (phase delay for single-frequency), so its variation under a parameter shift is

\delta U_n = P_n\,\delta\tau_n = \frac{P_n\,\delta\phi_n}{\omega}. \tag{4}

This identification — energy change equals power times phase-delay change — is the classical-wave analogue of the Wigner–Smith time-delay relation, here applied to a spatial parameter rather than frequency.

2.4 Connecting Port Quantities to the S-matrix

From the polar decomposition of $$b_n$$ ,

b_n^*\,\delta b_n = \tfrac{1}{2}\delta P_n + i\,P_n\,\delta\phi_n. \tag{5}

Summing over all ports and using matrix notation ( $\mathbf{b} = \mathbf{S}\mathbf{a}$ ):

\sum_n b_n^*\,\delta b_n = \mathbf{a}^\dagger\mathbf{S}^\dagger\delta\mathbf{S}\,\mathbf{a} = \frac{1}{2}\sum_n\delta P_n + i\sum_n P_n\,\delta\phi_n. \tag{6}

Lossless condition: Since $\mathbf{S}$ is unitary, total output power equals total input power. Under a virtual displacement (the incident amplitudes $\mathbf{a}$ are held fixed), $\sum_n \delta P_n = 0$ . The real part of (6) vanishes, leaving

\mathbf{a}^\dagger\mathbf{S}^\dagger\delta\mathbf{S}\,\mathbf{a} = i\sum_n P_n\,\delta\phi_n = i\,\omega\,\delta U. \tag{7}

2.5 The Central Result

Substituting (7) into the virtual work principle (2)–(3):

\sum_i \mathbf{F}_i\cdot\delta\mathbf{r}_i = -\delta U = \frac{1}{\omega}\,i\,\mathbf{a}^\dagger\mathbf{S}^\dagger\delta\mathbf{S}\,\mathbf{a} = \mathbf{a}^\dagger\mathbf{Q}_i\,\mathbf{a}\cdot\delta\mathbf{r}_i, \tag{8}

where the last equality uses definition (1) and the fact that $\mathbf{Q}_i$ is Hermitian (so its expectation value is real) in the lossless case.

Since the virtual displacements $\{\delta\mathbf{r}_i\}$ of distinct objects correspond to independent generalized coordinates (no geometric coupling between separate objects), they may be varied independently. Therefore, (8) holds component-by-component:

\boxed{\mathbf{F}_i = \langle \mathbf{Q}_i \rangle \equiv \mathbf{a}^\dagger \mathbf{Q}_i\,\mathbf{a}.} \tag{9}

This is the central equation of the GWS radiation-force framework, derived purely from virtual work and port-phase bookkeeping.

2.6 Remark on Independence and Multi-Object Optimization

The independence of $\{\delta\mathbf{r}_i\}$ in (8)–(9) has an immediate consequence for multi-objective optimization: there is no fundamental trade-off between simultaneously maximizing $\mathbf{F}_i$ and $\mathbf{F}_j$ for distinct objects $i \neq j$ arising from the virtual work structure itself. Any apparent trade-off is entirely a consequence of the specific scattering geometry (encoded in $\mathbf{S}$ ), not an intrinsic uncertainty principle. In particular:

If the system is block-diagonal (objects uncoupled), forces are simultaneously maximizable by construction.
Even in a coupled system, if $\mathbf{Q}_i$ and $\mathbf{Q}_j$ share a common maximum eigenvector, both forces are simultaneously maximized — irrespective of whether $[\mathbf{Q}_i, \mathbf{Q}_j] = 0$ or not.

A Robertson-type bound on $\langle\mathbf{Q}_i\rangle\langle\mathbf{Q}_j\rangle$ constrains variances, not maxima, and therefore does not constitute an uncertainty principle for radiation-force optimization.

3. Extension to Lossy Scattering Systems

3.1 Modified S-matrix Structure

When the scattering medium is lossy, $\mathbf{S}$ is sub-unitary:

\mathbf{S}^\dagger\mathbf{S} = \mathbf{I} - \mathbf{A}, \qquad \mathbf{A} \geq 0, \tag{10}

where $\mathbf{A}$ is the positive semi-definite absorption matrix. Total output power is no longer conserved:

\sum_n \delta P_n = -\delta P_\mathrm{abs}, \tag{11}

where $\delta P_\mathrm{abs} > 0$ is the additional power absorbed when the object is displaced to $\mathbf{r}_i + \delta\mathbf{r}_i$ .

By recovering the DoFs of loss channels, i.e. the bath, we can have a better understanding of this fact. Now consider a full scattering matrix

\mathbf{S} = \begin{pmatrix} \mathbf{S}_s & \mathbf{C}\ \mathbf{D} & \mathbf{S}_b \end{pmatrix}

\begin{pmatrix} \mathbf{S}_s & \mathbf{C}\ \mathbf{D} & \mathbf{S}_b \end{pmatrix} \begin{pmatrix} \mathbf{S}_s^\dagger & \mathbf{D}^\dagger\ \mathbf{C}^\dagger & \mathbf{S}_b^\dagger \end{pmatrix} = \begin{pmatrix} I & 0\ 0 & I \end{pmatrix}

\mathbf{S}_s \mathbf{S}_s^\dagger + \mathbf{C}\mathbf{C}^\dagger = I

\mathbf{S}_s \mathbf{D}^\dagger + \mathbf{C}\mathbf{S}_b^\dagger = 0\Rightarrow \mathbf{C}\approx-\mathbf{S}_s\mathbf{D}^\dagger\mathbf{S}_b\Rightarrow \mathbf{C}\mathbf{C}^\dagger \approx \mathbf{S}_s\mathbf{D}^\dagger \mathbf{D}\mathbf{S}_s^\dagger

3.2 Force Decomposition

Returning to (6), the real part no longer vanishes. Defining the (now non-Hermitian) generalized GWS operator ( $\mathrm{G^2WS}$ 😂)

\tilde{\mathbf{Q}}_i = \,\mathbf{S}^\dagger\,\partial_i\mathbf{S}, \tag{12}

equation (6) gives

\mathbf{a}^\dagger\tilde{\mathbf{Q}}_i\,\mathbf{a} = -\frac{1}{2}\frac{\partial P_\mathrm{abs}}{\partial \mathbf{r}_i} + i\,\omega\frac{\partial U}{\partial \mathbf{r}_i}. \tag{13}

The virtual work should also be corrected

\delta W + \delta W_\text{abs} = -\delta U

F_i \cdot \delta r_i +\delta W_\text{abs} = -\sum_n P_n \delta\phi_n /\omega= a^\dagger Q_i a/i\omega+\frac{1}{2i\omega}\frac{\partial P_\text{abs}}{\partial r_i}

which then yields a two-component force:

\mathbf{F}_i = \underbrace{\mathrm{Im}\!\left(\langle\tilde{\mathbf{Q}}_i\rangle\right)}_{\text{reactive (radiation pressure)}} + \underbrace{\frac{1}{2\omega}\,\frac{\partial P_\mathrm{abs}}{\partial \mathbf{r}_i}}_{\text{dissipative (absorption force)}}. \tag{15}

Decomposing $\tilde{\mathbf{Q}}_i = \mathbf{H}_i + i\mathbf{K}_i$ into Hermitian parts ( $\mathbf{H}_i^\dagger = \mathbf{H}_i, \mathbf{K}_i^\dagger = \mathbf{K}_i$ ):

$\mathbf{H}_i = \tfrac{1}{2}(\tilde{\mathbf{Q}}_i + \tilde{\mathbf{Q}}_i^\dagger)$ encodes the absorption gradient: $\langle\mathbf{H}_i\rangle \propto -\partial_i P_\mathrm{abs}$ .
$\mathbf{K}_i = \tfrac{1}{2i}(\tilde{\mathbf{Q}}_i - \tilde{\mathbf{Q}}_i^\dagger)$ encodes the radiation force: $\mathbf{F}_i = \langle\mathbf{K}_i\rangle$ .

Note that $\mathbf{K}_i$ is Hermitian, so the radiation-force optimization problem

\max_{\mathbf{a}:\,\|\mathbf{a}\|=1}\;\langle\mathbf{K}_i\rangle \tag{16}

remains an eigenvalue problem — the maximum is the largest eigenvalue of $\mathbf{K}_i$ , achieved by its corresponding eigenvector. No nonlinear solver is required even in the lossy case.

3.3 Physical Interpretation

The two terms in (14) have natural counterparts in known physics:

Term	Electromagnetic	Acoustic
Reactive (Im part)	Radiation pressure, gradient force	Acoustic radiation force
Dissipative (Re part)	Photothermal / photophoretic force	Acoustic streaming drag

The virtual work framework unifies these mechanisms within a single port-based expression, whereas they are traditionally treated by separate theories (Maxwell stress tensor + heat equation; Gorkov potential + Navier–Stokes streaming). In the lossless limit $\mathbf{A} \to 0$ , the dissipative term vanishes and (14) reduces to (9).

3.4 Consequences for Multi-Objective Optimization in Lossy Media

In the lossy case, each objective function is $f_i = \langle\mathbf{K}_i\rangle$ , which remains a Hermitian quadratic form. The Pareto structure of the multi-objective problem is therefore governed by the joint spectral properties of $\{\mathbf{K}_i\}$ , not of the non-Hermitian $\{\tilde{\mathbf{Q}}_i\}$ . The remarks of Section 2.6 apply unchanged: simultaneous maximization is possible whenever the matrices $\{\mathbf{K}_i\}$ share a common maximum eigenvector, and no Robertson-type bound on the maxima follows from non-commutativity alone.

4. Relation to Smith (1960) and Wigner (1955)

The identification in equation (4), $\delta U_n = P_n\,\delta\phi_n/\omega$ , is the spatial-parameter analogue of Smith’s time-delay formula. Smith (1960) defined the lifetime matrix $\mathbf{Q}$ through the residence-time argument with $\partial_\omega \mathbf{S}$ ; the present derivation replaces $\partial_\omega$ by $\partial_{\mathbf{r}_i}$ and replaces the quantum mechanical flux normalization argument by the classical virtual work balance. The two derivations are structurally identical, confirming that the GWS radiation-force relation is the spatial-domain counterpart of the Wigner–Smith time-delay relation, both rooted in the same energy-phase bookkeeping.

This observation implies that the “central equation” of recent GWS-based radiation-force frameworks is not a new variational identity, but rather a direct consequence of Smith’s 1960 construction applied to spatial rather than spectral parameters — a fact made transparent by the virtual work approach but obscured in field-theoretic derivations.

5. Summary

Property	Lossless	Lossy
S-matrix	Unitary	Sub-unitary
GWS operator	Hermitian $\mathbf{Q}_i$	Non-Hermitian $\tilde{\mathbf{Q}}_i$
Radiation force	$\mathbf{F}_i = \langle\mathbf{Q}_i\rangle$	$\mathbf{F}_i = \langle\mathbf{K}_i\rangle$ (Im part)
Absorption force	—	$\propto \langle\mathbf{H}_i\rangle$ (Re part)
Optimization	Eigenvalue problem	Eigenvalue problem ( $\mathbf{K}_i$ )
Multi-object trade-off	Geometry-dependent, not intrinsic	Same

The virtual work derivation is shorter, more transparent, and more general than field-theoretic approaches. It makes the assumptions explicit (lossless $\leftrightarrow \sum\delta P_n = 0$ ), reveals the two-force decomposition in lossy media, and demonstrates that analytical eigenvalue methods suffice throughout — eliminating the need for general-purpose nonlinear solvers.

Note: Equations are written in the +j\omega t time-harmonic convention, consistent with the GWS literature.