Coupled Oscillators and Absorption

The spectrum of a coupled oscillator system

This note records another calculation in my quest to understand why water is blue. Last time we looked at absorption of light by a single oscillator; now we graduate to systems of coupled oscillators — systems closer to the stretching and bending vibrational modes in a molecule.

Again, we proceed in two stages: first the linear story which we can solve exactly, then the (simplest) nonlinear story, which we explore numerically.

Linear coupling

Linear coupling means a quadratic potential. Every quadratic form on $\mathbb R^n$ can be written

$$V(x) \;=\; \tfrac12\, x^T H x$$

for a unique real symmetric $H$. We’re studying vibrational modes, so $H$ is positive definite — the origin is a nondegenerate minimum.

The conservative equation of motion is

$$\ddot x \,+\, H x \;=\; 0.$$

A live view of $\ddot x + Hx = 0$ in the original $x$ coordinates. Edit any cell of the symmetric matrix $H$ above to reshape the system; the eigenvalues update on every keystroke. Click and drag any rail to set its initial position. Dotted gray lines connect rails $i$ and $j$ whenever $H_{ij} \ne 0$, with thickness scaling with $|H_{ij}|$ — energy sloshes between coupled rails as the modes beat against one another.

These systems can be solved exactly by decoupling. Orthogonally diagonalize as $H = O\Lambda O^T$ with $\Lambda = \mathrm{diag}(\omega_1^2, \ldots, \omega_n^2)$ (eigenvalues positive by positive-definiteness). In the eigenbasis $y = O^T x$ the equation decouples mode by mode:

$$\ddot y_i \,+\, \omega_i^2\, y_i \;=\; 0.$$

A coupled linear conservative system is $n$ independent simple harmonic oscillators in disguise.

The diagonalization made visible. Left: motion in $x$ coordinates, with dotted lines marking the nonzero couplings of $H$. Right: the same trajectory in modal coordinates $y = O^T x$ — $n$ independent simple harmonic oscillators each at $\omega_i$, no couplings. Push a single $x$-rail and watch energy distribute across multiple $y$-rails; the “Mode $i$” presets do the inverse, exciting one mode in $y$ and showing the corresponding pattern across multiple $x$-rails.

Driving and damping

What we’re really after isn’t the free theory but a system of oscillators driven by a time-dependent external force and damped by their environment:

$$\boxed{ \ddot x \,+\, 2\gamma\, \dot x \,+\, H x \;=\; A \cos(\Omega t), }$$

where $\gamma > 0$ is the (uniform) damping rate and $A \in \mathbb R^n$ the drive vector. The quantity we want is the absorbed power — the time-averaged dissipation rate $2\gamma \sum_i \dot x_i^2$ over one drive period $T = 2\pi/\Omega$,

$$P_{\text{abs}}(\Omega) \;=\; \frac{2\gamma}{T} \int_0^T \sum_{i=1}^n \dot x_i(t)^2 \, dt.$$

To compute this we first solve the system. The same diagonalization works: setting $y = O^T x$ and $\tilde A := O^T A$,

$$\ddot y \,+\, 2\gamma\, \dot y \,+\, \Lambda y \;=\; \tilde A \cos(\Omega t),$$

which is just $n$ independent driven damped harmonic oscillators

$$\ddot y_i \,+\, 2\gamma\, \dot y_i \,+\, \omega_i^2\, y_i \;=\; \tilde A_i \cos(\Omega t),$$

each solved in the last post and seen to be a transient decaying to a periodic steady state

$$y_i^{\text{ss}}(t) \;=\; \frac{\tilde A_i}{\sqrt{(\omega_i^2 - \Omega^2)^2 + 4\gamma^2\Omega^2}}\,\cos(\Omega t - \varphi_i).$$

(Decoupling depends on $\gamma$ being scalar. If different oscillators damp at different rates — say because they have different masses — the change of basis no longer decouples the system, and we are stuck with legitimately coupled modes. A different story.)

The absorption spectrum

Coming back to what we want: $P_{\text{abs}}(\Omega) = 2\gamma \langle |\dot x|^2 \rangle$, with $\langle \cdot \rangle = \tfrac1T \int_0^T \cdot \, dt$ the time average. Since $O$ is an isometry,

$$|\dot x|^2 \;=\; |O\dot y|^2 \;=\; |\dot y|^2,$$

so by linearity of $\langle \cdot \rangle$,

$$\begin{aligned} P_{\text{abs}}(\Omega) &\;=\; 2\gamma \langle |\dot y|^2 \rangle \\ &\;=\; 2\gamma \left\langle \sum_{i=1}^n \dot y_i^2 \right\rangle \\ &\;=\; 2\gamma \sum_{i=1}^n \langle \dot y_i^2 \rangle \\ &\;=\; \sum_{i=1}^n P_{\omega_i}(\Omega), \end{aligned}$$

where $P_{\omega_i}(\Omega) := 2\gamma \langle \dot y_i^2 \rangle$ is the per-mode absorption from the previous post. Plugging in the steady-state response $y_i^{\text{ss}}$ gives the explicit form,

$$\boxed{ P_{\text{abs}}(\Omega) \;=\; \sum_{i=1}^n \frac{\gamma\, \tilde A_i^{\,2}\, \Omega^2}{(\omega_i^2 - \Omega^2)^2 + 4\gamma^2 \Omega^2}. }$$

The absorption spectrum as a sum of single-mode contributions. The bold curve is the total $P_{\text{abs}}(\Omega)$; the faint curves underneath are the individual $P_{\omega_i}$. Drag $\gamma$ to broaden or narrow every peak together; edit $H$ to slide their positions and adjust their relative weights $\tilde A_i^{\,2} = (O^T A)_i^{\,2}$. The drive $A = (1, 1, \ldots, 1)^T$ is uniform across all oscillators in this demo.

The linear story is the previous post, $n$ times over: a sum of single-oscillator spectra, weighted by the drive’s projection onto each mode.

Nonlinear coupling

Our next goal is to understand nonlinear coupled systems, where we can’t solve things exactly and must resort to numerics. As in the previous post, we’ll do this by studying the simplest possible nonlinear system.

We work with two oscillators $x_1, x_2$ at natural frequencies $\omega_1, \omega_2$ and look for the simplest cross-coupling to add to the harmonic potential $\tfrac12 \omega_1^2 x_1^2 + \tfrac12 \omega_2^2 x_2^2$. Quadratic cross-terms are absorbed by the linear theory, so we need cubic. Of the four cubic monomials $x_1^3, x_2^3, x_1^2 x_2, x_1 x_2^2$, the pure cubes don’t couple, and the two genuine cross-terms are equivalent under $x_1 \leftrightarrow x_2$. Without loss of generality

$$V(x_1, x_2) \;=\; \tfrac12\, \omega_1^2 x_1^2 \,+\, \tfrac12\, \omega_2^2 x_2^2 \,+\, \epsilon\, x_1^2 x_2.$$

This $V$ is unbounded below — same caveat as the cubic case last post — with a saddle at finite distance opening an escape route in $x_1$ once $x_2 < -\omega_1^2/(2\epsilon)$.

With uniform damping $\gamma$ and a common drive $A\cos(\Omega t)$ on each oscillator, the equations of motion are

$$\boxed{ \begin{aligned} \ddot x_1 \,+\, 2\gamma\, \dot x_1 \,+\, \omega_1^2\, x_1 \,+\, 2\epsilon\, x_1 x_2 &\;=\; A\cos(\Omega t), \\ \ddot x_2 \,+\, 2\gamma\, \dot x_2 \,+\, \omega_2^2\, x_2 \,+\, \epsilon\, x_1^2 &\;=\; A\cos(\Omega t). \end{aligned} }$$

The coupling is asymmetric: $\ddot x_1$ has its stiffness modulated by $x_2$, while $\ddot x_2$ feels a force quadratic in $x_1$.

The system in motion. Each $x_i$ rides its own rail, connected by a line marking the coupling $\epsilon x_1^2 x_2$. Slate-blue driver balls visualize the common drive $A\cos(\Omega t)$ on each rail. Click and drag to set initial conditions; sliders control $\gamma, A, \Omega, \omega_1, \omega_2, \epsilon$. After a transient, trajectories typically settle to periodic motion locked to the drive — though at large $A$ or $\epsilon$ they can escape over the saddle.

The same dynamics in phase space show clearly the existence of limiting periodic behavior, after a transient initial component dies away. The demo below shows both, in the $(x_1, x_2)$ plane:

Light gray traces the integrated trajectory from the user-chosen IC; maroon overlays the eventual periodic limit cycle, detected automatically. Click anywhere to set the IC; sliders update live.

Finding the maroon curve is a Poincaré check. Linear transients decay as $e^{-\gamma t}$, so integrating from rest for $t \gg 1/\gamma$ puts us comfortably on the steady state. We then test whether the state has returned at successive integer multiples of the drive period $T = 2\pi/\Omega$; the smallest such integer $k$ is the period of the limit cycle, and integrating $kT$ forward traces the closed curve.

The absorption spectrum

The same integral defines absorption,

$$P_{\text{abs}}(\Omega) \;=\; \frac{2\gamma}{T} \int_0^T \bigl(\dot x_1^2 + \dot x_2^2\bigr) \, dt,$$

now evaluated on the numerically-computed steady state. For each $\Omega$ we integrate from rest for $M$ drive periods chosen so $\gamma M T \geq 40$ — the residual transient $\sim e^{-\gamma M T}$ is then well below numerical noise — and time-average $2\gamma(\dot x_1^2 + \dot x_2^2)$ over a few drive periods to get $P_{\text{abs}}(\Omega)$.

The absorption spectrum on a log-$y$ scale. The dashed curve is the analytic $\epsilon = 0$ baseline; the solid curve is the full nonlinear $P_{\text{abs}}$, built up point by point as the sweep progresses. Vertical guides mark $\omega_1$ and $\omega_2$; gray bands flag drive frequencies where the system escapes and no steady state exists. Sliders control $\gamma, \epsilon, \omega_1, \omega_2$.

A few features stand out. The two main peaks at $\omega_1$ and $\omega_2$ — the linear answer — survive, slightly distorted by the nonlinearity. Below them, smaller peaks appear as $\epsilon$ grows, sitting roughly at $\omega_i/n$ for small integers $n$ — the same kind of subharmonic structure we saw in the previous post’s cubic case. The gray bands flag where the unbounded well kicks in: at large $\epsilon$ or $A$ some trajectories escape and there is no steady state.

Real spectra are richer still: water’s, for instance, has prominent overtones at integer multiples $n\omega_i$ above the fundamentals, and combination bands at sum frequencies like $\omega_1 + \omega_2$. Neither appears in the picture above, which shows only the fundamentals and a subharmonic ladder below them. I would guess this means we need to think harder about what is really going on, before we are ready to model water.