With the operator \mathbf{A}(\omega) of the perturbation network assembled, every acoustic question becomes a statement about that one matrix. Four analyses read different things from it: its resonances and their growth rates, an open-loop stability count taken on the real-frequency axis, the flux of acoustic energy and what it reveals about passivity, and the field forced by a prescribed excitation.
1.1 The sign convention
The acoustics are posed with the time dependence X'(t) = \Re\{\widehat{X}\,e^{\mathrm{i}\omega t}\}, and a free mode of the network has a complex frequency \omega = \omega_r + \mathrm{i}\omega_i. Its time evolution is then e^{\mathrm{i}\omega t} = e^{\mathrm{i}\omega_r t}\,e^{-\omega_i t}, so a positive imaginary part decays and a negative one grows. The modal growth rate is therefore defined as:
\sigma \;=\; -\omega_i \;=\; -\Im(\omega),
\qquad\text{a mode is unstable when } \sigma > 0 \ (\Im(\omega) < 0),
where \omega_r/(2\pi) is the modal frequency in hertz and \sigma the growth rate in inverse seconds. It should be emphasized that this convention — growth as the negative imaginary part — is the authoritative one throughout the implementation, fixed by a lossy-duct test that pins a passive decaying mode to \Im(\omega) > 0; it is stated here once and used without further comment below.
1.2 Modal stability
The resonances of the network are the frequencies at which the operator is singular, so the modal problem is the nonlinear eigenproblem, given as:
where each root \omega yields a modal frequency \omega_r/(2\pi) and a growth rate \sigma = -\omega_i. The operator is entire in \omega — its only \omega-dependence is through \mathrm{i}\omega\mathbf{M}, the duct phases e^{-\mathrm{i}\omega\tau}, and any analytic source or boundary transfer function, none of which has a pole — so the roots are isolated and can be counted.
The roots are found by the contour-integral method of Beyn (Beyn 2012), which recovers all eigenvalues inside a chosen contour from moments of the resolvent, given as:
where \widehat{\mathbf{V}} is a random probe block and the contour is an ellipse spanning the frequency band and growth range of interest. The rank of the zeroth moment \mathbf{A}_0 is the number of enclosed modes, and a small dense eigenproblem built from the two moments returns the eigenvalues and eigenvectors — the method never forms \det\mathbf{A} itself, which would overflow.
That rank cannot be read off \mathbf{A}_0 alone. When the contour encloses no eigenvalue the integrand is analytic, \mathbf{A}_0 vanishes identically, and its computed singular values are nothing but quadrature error; a threshold placed relative to the largest of them then reports the full probe width where the true rank is zero, and the search returns as many modes as it was given probes. Nor is there an absolute threshold that separates the two cases, because the floor is set by the truncation error from poles lying just outside the contour rather than by round-off.
The rank is therefore supplied from outside the method. The argument principle counts the enclosed roots by the winding of the determinant, given as:
N \;=\; \frac{1}{2\pi\mathrm{i}}\oint \frac{\det'\mathbf{A}(z)}{\det\mathbf{A}(z)}\,\mathrm{d}z \;=\; \frac{1}{2\pi}\,\Delta\arg\det\mathbf{A},
where the phase of \det\mathbf{A} is accumulated from the triangular factors of an LU decomposition so that its magnitude is never formed. The search region is covered by overlapping elliptical sub-contours; on each, N is evaluated first and handed to Beyn as the rank of its moment matrix, so a sub-contour enclosing nothing is skipped rather than mined for noise. The same count over the region as a whole then certifies the mode set as complete rather than merely “some modes found”, and the driver re-tiles and re-searches until the two agree. The sub-contours must genuinely cover the region they certify: ellipses laid side by side pinch at their seams, and a strongly damped mode sitting in such a gap would be counted yet never found.
The winding is trustworthy only where the counting contour resolves it: the phase is accumulated increment by increment, and each per-step rotation of \arg\det\mathbf{A} must remain below \pi for it to unwrap onto the correct branch. A spectrum denser than the contour can resolve folds the phase by more than \pi between neighbouring nodes, and the winding then aliases onto a spurious integer cleanly enough that its distance to the nearest integer betrays nothing; only the largest per-step rotation, as it approaches \pi, reveals the aliasing. This is the regime of a low-Mach flowing network, whose convected entropy and composition modes are spaced by the long transit time \tau_u = L/\overline{u} rather than the acoustic \tau_\pm = L/(\overline{c}\pm\overline{u}), and so crowd the region far more densely than the acoustic count that sizes the contour anticipates. There the completeness certificate is withheld rather than asserted from an aliased count, and the acoustic spectrum is recovered instead under the isentropic reduction, which pins the convected wave to zero (tests: test_low_mach_dense_spectrum_is_not_falsely_certified, test_low_mach_choked_nozzle_warns_not_crashes).
Membership must also be decided on a scale-invariant footing. A network operator assembles rows in incompatible units — pressure, velocity, mass flow — so \operatorname{cond}\mathbf{A} routinely reaches 10^{12} at every frequency, and a residual normalised by \max|\mathbf{A}| is satisfied by a near-null vector wherever one cares to look. Modes are accepted on the residual of the equilibrated operator \mathbf{A}_s = \mathbf{D}_r\mathbf{A}\mathbf{D}_c, whose diagonal scalings are frozen once at the band centre; the Newton polish that precedes the test converges on |\Delta\omega|, which needs no scale at all.
Each mode carries an eigenvector, projected per edge to the characteristic amplitudes (f, g, h) through \mathbf{L}_e (see characteristics), which is the mode’s spatial wave pattern. The behaviour is illustrated in Figure 1 for a Rijke tube — a duct with a compact heat source — whose modes move across the stability boundary as the flame is activated (tests: test_beyn_oracle_lossy_complex_modes_and_sign, test_eigenmodes_certified_count_matches, test_beyn_moment_rank_is_ambiguous_on_an_empty_contour, test_no_modes_survive_an_eigenvalue_free_region, test_n_tau_flame_drives_self_excited_instability).
The frequency and growth bands set the semi-axes of the search ellipse, not the sides of a rectangle: a mode near a corner of the implied box lies outside the region, and is neither counted nor returned.
Figure 1: Eigenmode spectrum of a Rijke tube — growth rate against modal frequency — for three flame settings, each a certified contour-eigensolver search of \det\mathbf{A}(\omega)=0. Passive (no unsteady heat release): the modes are the tube’s acoustic resonances, essentially neutral. Stabilizing flame (n=0.8,\ \tau=1.5\,ms): the fundamental is pushed further below the boundary. Destabilizing flame (n=0.8,\ \tau=4\,ms): a mode is lifted above the growth-rate boundary (dashed) into instability. Points above the line grow; points below decay.
The contour eigensolver is powerful but not universal, and three circumstances defeat it: a measured flame transfer function is defined only on the real axis and cannot be evaluated on the complex contour; a flowing network fills its spectrum with dense, near-marginal convected entropy and composition modes that overwhelm the completeness certificate; and a long convective delay makes e^{-\mathrm{i}\omega\tau_u} overflow once the contour reaches into the complex plane. For these regimes stability is assessed on the real axis, by a Nyquist argument (Nyquist 1932) on the feedback the source introduces.
The construction splits the operator into its passive part and the source, given as:
where \mathbf{A}_0 is the network with every source switched off and \mathbf{S} is the low-rank source, one rank-one term per feedback with transfer function \mathcal{F}_k, injection vector \mathbf{a}_k, and sensing vector \mathbf{b}_k. The matrix-determinant lemma then factors the determinant into the passive one and a small determinant of the return ratio, given as:
where \mathbf{L} is the return ratio (a scalar for a single flame) and D(\omega) = \det(\mathbf{I} - \mathbf{L}) = \det\mathbf{A}/\det\mathbf{A}_0 is the stability determinant. Evaluated along the real-frequency axis, the locus of D(\omega) encircles the origin once for each unstable mode, so the count follows from the encirclements, given as:
n_{\text{unstable}} \;=\; -\tfrac{1}{2}\,(\text{winding of } D \text{ about the origin}),
where the factor of one-half accounts for the closed contour including the negative-frequency image D(-\omega) = \overline{D(\omega)}, and the least value \min_\omega|D| is the stability margin, vanishing at the onset of instability. An important qualification is that this counts unstable modes relative to\mathbf{A}_0, so it equals the absolute count only when the passive network is itself stable, a condition the driver checks with a rational fit of D, and the count is limited to the maximum swept frequency (tests: test_unstable_count_matches_eigenmodes, test_entropy_path_is_the_sole_destabilizer, test_tabulated_ftf_supported_on_the_real_axis).
1.3 Acoustic energy and passivity
A complementary reading of the operator is energetic, through the mean-flow acoustic energy of Myers (Myers 1991). For a plane-wave field the downstream energy flux (intensity) and the energy density are given as:
I = \tfrac{1}{2}\overline{\varrho}\,\overline{c}\Big[(1 + M)^2|\widehat{f}|^2 - (1 - M)^2|\widehat{g}|^2\Big],
\qquad
e = \tfrac{1}{2}\overline{\varrho}\Big[(1 + M)|\widehat{f}|^2 + (1 - M)|\widehat{g}|^2\Big],
where \widehat{f} and \widehat{g} are the downstream and upstream acoustic amplitudes and M the mean Mach number, and for a single downstream wave I/e = \overline{u} + \overline{c}, the transport of energy at the group speed. These give a sharp statement of passivity: an energy-neutral termination reflects with a magnitude fixed by the mean flow, given as:
and a reflection magnitude above this bound adds acoustic energy — an active boundary that behaves as a source. The energy view also furnishes an independent check on the eigensolver’s growth rates: for a source-free mode the stored energy obeys \mathrm{d}E/\mathrm{d}t = 2\sigma E, which equals the net acoustic power delivered through the boundaries, so the sign of the net boundary power must match the sign of the growth rate — a physically grounded cross-check the code performs mode by mode (tests: test_passive_bound_is_the_zero_flux_reflection, test_boundary_power_sign_matches_growth_every_mode, test_modal_energy_balance_recovers_growth_rate).
1.4 Forced response
The final analysis drives the network rather than letting it oscillate freely. With a forcing prescribed at a driven terminal, the response is one sparse solve per frequency, given as:
where \widehat{\mathbf{b}} is the forcing vector, nonzero only at the driven terminal rows. From the solved nodal field one reads the per-edge wave amplitudes, the reflection \widehat{g}/\widehat{f} at any station, and the stored acoustic energy. Because the system is linear, the response scales exactly with the forcing amplitude and superposes over several drives, and it becomes singular precisely at the resonances where \det\mathbf{A}(\omega) = 0 — the same modes the eigensolver locates (tests: test_forced_response_scales_linearly_with_amplitude, test_forced_response_obeys_superposition).
All four analyses read the single operator \mathbf{A}(\omega), and all four presume that every element in it is modeled. When one element is not, e.g. a transfer matrix or a transfer function is unknown, that dynamic response can be recovered from a measured network response by de-embedding, which is identification.