7.7 Applications to springs, circuits, epidemics and qubits

Four real systems run on the machinery of Unit 7. Coupled masses split into independent oscillators along the eigenvectors of a symmetric matrix. An RLC circuit is a damped oscillator, so the same equation, the same roots and the same resonance. The SIR epidemic has a threshold R0 that is the sign of one eigenvalue at the disease-free equilibrium. The state of a two-level quantum system evolves by the matrix exponential of -iH/hbar.

Unit 7 produced tools: solve $y'' + by' + cy = R$ (7.2, 7.3), solve $x' = Ax$ by the matrix exponential (7.4), and read stability off eigenvalues or off the Jacobian (7.6). This chapter points them at four real systems, each with numbers. The pattern is the same every time. A law of nature (Newton, Kirchhoff, “infections need a meeting between an infected and a susceptible person”, Schrödinger) gives a system of differential equations. The system is linear, or becomes linear near an equilibrium. And then the eigenvalues carry the meaning: imaginary parts are frequencies, real parts are decay and growth rates, and the sign of one eigenvalue can be the difference between an epidemic and none.

Before you start

You need the spectral theorem for symmetric matrices from 3.1, the roots of the characteristic polynomial from 7.2, the $2\times2$ matrix exponential from 7.4, and the language of equilibria and stability from 7.6. Three quick checks:

Quick checkFor a real symmetric matrix, eigenvectors that belong to different eigenvalues are

Quick checkThe equation

x'' + 2x' + 5x = 0

describes

Quick checkIf

A^2 = -\omega^2I

for a

2 \times 2

matrix

A

and a number

\omega > 0

, what is

e^{tA}

The picture

Two identical masses hang between walls and are joined to each other by a weaker spring. Press Play to watch the time slider run. Start with Push mass 1 only: the first mass swings, slows down as the second one picks up, and after about five seconds the first mass is almost still while the second swings wide; then the energy flows back. Now choose Mode 1 or Mode 2. Each of these starts both masses in a special pattern, and then both masses oscillate at a single frequency for ever, without any exchange of energy. These special patterns are the normal modes, and every motion is a sum of them.

Where this came from· Kermack and McKendrick, Rabi and Schrödinger

In 1926 Schrödinger published his wave equation in a series of four papers titled Quantisierung als Eigenwertproblem in the Annalen der Physik: quantisation as an eigenvalue problem, the energies being eigenvalues of an operator. In 1927 W. O. Kermack and A. G. McKendrick published “A contribution to the mathematical theory of epidemics” in the Proceedings of the Royal Society A (volume 115, pages 700 to 721). Their paper established a threshold: for given rates of infection and recovery, a single infected person introduced into a population starts an epidemic only if the density of susceptible people is above a critical value. In 1937 Isidor Rabi published “Space quantization in a gyrating magnetic field” in Physical Review (volume 51, page 652), computing how a spin responds to a rotating magnetic field. That calculation became the theoretical basis for magnetic resonance, from nuclear magnetic resonance to the MRI scanner.

The idea, in words

Why· Find the eigen-directions, and each part of the problem becomes a single number

All four systems are, after a change of variables, copies of the scalar equations you know. Coupled masses: the symmetric matrix has an orthonormal basis of eigenvectors, and in that basis the system falls apart into independent oscillators $c'' = -\omega^2c$ . The circuit: a single oscillator with damping, so a complex pair of roots, the same equation as a mass on a spring. The epidemic: a nonlinear system, linearized at the state where nobody is infected, where one eigenvalue is $\gamma(R_0 - 1)$ and its sign is the threshold. The qubit: a matrix exponential whose matrix has imaginary eigenvalues, so the state rotates and never decays. In every case the question ‘what does the system do?’ becomes ‘what are the eigenvalues, and what do they mean?’

Definitions and theorems

Coupled springs and normal modes

Take $n$ masses on a line, each joined by springs to its neighbours or to the walls. Let $x_i$ be the displacement of mass $i$ from its resting place. Hooke’s law says a spring pulls with force $k\times(\text{extension})$ , and Newton’s law says $m x_i''$ is the total force. For two masses $m$ between walls, with a spring $k$ from each wall and a coupling spring $k_c$ between them,

m\,x_1'' = -k\,x_1 + k_c(x_2 - x_1), \qquad m\,x_2'' = -k\,x_2 - k_c(x_2 - x_1).

In matrix form this is $x'' = -Mx$ with

M = \frac1m\begin{bmatrix}k + k_c & -k_c\\-k_c & k + k_c\end{bmatrix}.

The matrix $M$ is symmetric. It is also positive definite: the stored energy $\frac12kx_1^2 + \frac12k_c(x_2 - x_1)^2 + \frac12kx_2^2$ is positive for every $x \ne 0$ , and this energy is $\frac m2x^{\mathsf T}Mx$ (3.2). The same shape $x'' = -Mx$ appears for any number of masses with equal masses $m$ .

Definition B

Normal mode

A normal mode of $x'' = -Mx$ is a solution of the form $x(t) = \cos(\omega t - \phi)\,v$ with a fixed non-zero vector $v$ (the mode shape) and numbers $\omega > 0$ (the natural frequency) and $\phi$ . Substituting, $-\omega^2\cos(\omega t - \phi)v = -\cos(\omega t - \phi)Mv$ , so $Mv = \omega^2v$ : the mode shapes are the eigenvectors of $M$ and the natural frequencies are the square roots of its eigenvalues. For $M$ positive definite all eigenvalues are positive, so all frequencies are real.

In one lineA solution of $x'' = -Mx$ in which every coordinate oscillates at the same frequency in a fixed pattern: $x(t) = \cos(\omega t - \phi)\,v$ . Then $v$ is an eigenvector of $M$ and $\omega^2$ its eigenvalue.

Why it existsComplicated motion of many coupled parts is a sum of independent simple motions, one for each eigenvector. The frequencies of a structure, a molecule or an instrument are the square roots of the eigenvalues.

Theorem B standard; spectral theorem (3.1) with the solutions of 7.2

Normal-mode decomposition

Let $M$ be a real symmetric positive definite $n \times n$ matrix with an orthonormal basis of eigenvectors $v_j$ , $Mv_j = \omega_j^2v_j$ , $\omega_j > 0$ . Then the unique solution of $x'' = -Mx$ with given $x(0)$ and $x'(0)$ is

x(t) = \sum_{j=1}^{n}\Big(a_j\cos\omega_jt + b_j\sin\omega_jt\Big)v_j, \qquad a_j = v_j\cdot x(0),\ \ b_j = \frac{v_j\cdot x'(0)}{\omega_j}.

Each term is a normal mode, and the energy $\frac12\lVert x'\rVert^2 + \frac12x^{\mathsf T}Mx$ is the sum of the energies of the modes, each constant in time.

In one lineIf $M$ is symmetric positive definite with orthonormal eigenvectors $v_1, \dots, v_n$ and eigenvalues $\omega_j^2$ , then every solution of $x'' = -Mx$ is $x(t) = \sum_j\big(a_j\cos\omega_jt + b_j\sin\omega_jt\big)v_j$ with $a_j = v_j\cdot x(0)$ and $b_j = v_j\cdot x'(0)/\omega_j$ .

Why it is trueIn the basis $v_1, \dots, v_n$ the coordinates $c_j = v_j\cdot x$ satisfy $c_j'' = -\omega_j^2c_j$ , because $M$ is symmetric. Each is an ordinary oscillator.

Proofstandard

The idea: Expand x in the eigenvector basis. Use symmetry of M to move M across the dot product, and the coordinates decouple into n copies of the scalar oscillator of 7.2.

Write

x(t) = \sum_jc_j(t)v_j

with

c_j(t) = v_j\cdot x(t)

The spectral theorem (3.1) gives an orthonormal eigenbasis, and coordinates in an orthonormal basis are dot products (2.9).

c_j'' = -v_j\cdot Mx

Differentiate twice:

c_j'' = v_j\cdot x''

, since

v_j

is constant. Substitute

x'' = -Mx

c_j'' = -(Mv_j)\cdot x = -\omega_j^2\,v_j\cdot x = -\omega_j^2c_j

M

is symmetric, so

v_j\cdot Mx = (Mv_j)\cdot x

. And

Mv_j = \omega_j^2v_j

c_j(t) = a_j\cos\omega_jt + b_j\sin\omega_jt

The equation

c'' + \omega^2c = 0

has characteristic roots

\pm i\omega

, so its solutions are

a\cos\omega t + b\sin\omega t

(7.2). Here

a = c_j(0) = v_j\cdot x(0)

and

b\omega = c_j'(0) = v_j\cdot x'(0)

Hence the formula for

x(t)

Uniqueness of the solution of the initial value problem for the linear system (7.1): the formula is a solution with the right initial data, so it is the solution.

The total energy is a sum of constants, one per mode.

The energy:

\lVert x'\rVert^2 = \sum c_j'^2

and

x^{\mathsf T}Mx = \sum\omega_j^2c_j^2

(orthonormal basis). Each

\frac12c_j'^2 + \frac12\omega_j^2c_j^2 = \frac12\omega_j^2(a_j^2 + b_j^2)

is constant.

\blacksquare

Proof workshopRebuild: coupled oscillators split into normal modesstandard

Tap the lines in the order a proof would write them. Each line can only use what is already above it.

Worked example new numbers; compare 7.2 (y'' + ω²y = 0)

Two masses m = 0.5 kg, springs k = 20 N/m to the walls and kc = 2 N/m between them. Find the normal modes, then solve the motion when mass 1 is pulled 1 cm and released.

Solution.

Write $M$ . $M = \frac{1}{0.5}\begin{bmatrix}22 & -2\\-2 & 22\end{bmatrix} = \begin{bmatrix}44 & -4\\-4 & 44\end{bmatrix}$ (in $\text{s}^{-2}$ ). $k + k_c = 22$ , divide by $m = 0.5$ .
Eigenvalues and eigenvectors. $M\begin{bmatrix}1\\1\end{bmatrix} = \begin{bmatrix}40\\40\end{bmatrix}$ and $M\begin{bmatrix}1\\-1\end{bmatrix} = \begin{bmatrix}48\\-48\end{bmatrix}$ . So $\omega_1^2 = 40$ , mode shape $(1, 1)$ ; $\omega_2^2 = 48$ , mode shape $(1, -1)$ . Trace is $88 = 40 + 48$ and determinant $44^2 - 16 = 1920 = 40\cdot48$ .
Frequencies. $\omega_1 = \sqrt{40} = 6.3246$ rad/s ( $1.007$ Hz) and $\omega_2 = \sqrt{48} = 6.9282$ rad/s ( $1.103$ Hz). In mode 1 the masses move in step and the coupling spring is never stretched, so the frequency is the single-mass value $\sqrt{k/m} = 6.3246$ . In mode 2 they move in opposition and the coupling spring is stretched twice as much, giving $\sqrt{(k + 2k_c)/m}$ .
Starting data. $x(0) = (0.01, 0)$ m, $x'(0) = 0$ . The orthonormal modes are $\frac{1}{\sqrt2}(1, 1)$ and $\frac1{\sqrt2}(1, -1)$ , so $a_1 = a_2 = \frac{0.01}{\sqrt2}$ , $b_1 = b_2 = 0$ . $a_j = v_j\cdot x(0)$ and the masses start at rest.
The motion. $x_1 = \frac{0.01}{2}\big(\cos\omega_1t + \cos\omega_2t\big)$ and $x_2 = \frac{0.01}{2}\big(\cos\omega_1t - \cos\omega_2t\big)$ . $x = a_1\cos\omega_1t\cdot v_1 + a_2\cos\omega_2t\cdot v_2$ with $v_{1,2} = \frac1{\sqrt2}(1, \pm1)$ .
Recognise beats. With $\bar\omega = \frac{\omega_1 + \omega_2}{2} = 6.6264$ and $\Omega = \frac{\omega_2 - \omega_1}{2} = 0.3018$ , the sum-to-product identities give $x_1 = 0.01\cos(\bar\omega t)\cos(\Omega t)$ and $x_2 = 0.01\sin(\bar\omega t)\sin(\Omega t)$ . $\cos A + \cos B = 2\cos\frac{A+B}2\cos\frac{A-B}2$ and $\cos A - \cos B = -2\sin\frac{A+B}2\sin\frac{A-B}2$ (note the sign: $x_2$ with $\sin\bar\omega t\sin\Omega t$ ).
Read off the energy transfer. Mass 1 oscillates at about $6.63$ rad/s with an envelope $\cos(0.3018t)$ ; mass 2 with the envelope $\sin(0.3018t)$ . The envelope of mass 1 first vanishes when $0.3018t = \frac\pi2$ , that is $t = 5.20$ s: all the energy has moved to mass 2. The whole cycle takes $\frac{2\pi}{\omega_2 - \omega_1} = 10.4$ s. Neither frequency in the beating motion is a natural frequency of the system. The beat comes from the two natural frequencies being close, which happens when the coupling $k_c$ is weak.

Answer. Normal modes $(1, 1)$ at $6.3246$ rad/s and $(1, -1)$ at $6.9282$ rad/s; pulling one mass gives beats that transfer all the energy to the other mass every $5.20$ s.

Your answerIn the example above, after how many seconds has all the energy first passed from mass 1 to mass 2? (The time is

\pi/(\omega_2 - \omega_1)

with

\omega_1 = \sqrt{40}

\omega_2 = \sqrt{48}

Quick checkYou push only the first of two identical, weakly coupled oscillators and let go. What happens?

Worked example standard; new numbers

Three identical masses in a row between two walls (all springs k, all masses m). Find the normal modes and their frequencies.

Solution.

The matrix. $M = \frac km\begin{bmatrix}2 & -1 & 0\\-1 & 2 & -1\\0 & -1 & 2\end{bmatrix}$ . Mass $i$ feels $-k(x_i - x_{i-1}) - k(x_i - x_{i+1})$ : $2$ on the diagonal, $-1$ for each neighbour; the walls hold $x_0 = x_4 = 0$ .
Test the vector $(1, 0, -1)$ . $M(1, 0, -1)^{\mathsf T} = \frac km(2, 0, -2)^{\mathsf T} = 2\frac km(1, 0, -1)^{\mathsf T}$ . Eigenvalue $2k/m$ . Row 2 is $-1\cdot1 + 2\cdot0 - 1\cdot(-1) = 0$ .
Test $(1, \pm\sqrt2, 1)$ . $M(1, s, 1)^{\mathsf T} = \frac km(2 - s,\ -2 + 2s,\ 2 - s)^{\mathsf T}$ , and this equals $\lambda(1, s, 1)$ when $\lambda = (2 - s)\frac km$ and $\lambda s = (2s - 2)\frac km$ , that is $s(2 - s) = 2s - 2$ , so $s^2 = 2$ . Hence $s = \pm\sqrt2$ with eigenvalues $(2 \mp \sqrt2)\frac km$ . Require the same $\lambda$ from rows 1 and 2.
Frequencies. $\omega^2 = (2 - \sqrt2)\frac km$ , $2\frac km$ , $(2 + \sqrt2)\frac km$ . With $k = 20$ N/m and $m = 0.5$ kg ( $k/m = 40$ ): $\omega = 4.841,\ 8.944,\ 11.686$ rad/s. $\sqrt{40\times0.5858} = 4.841$ , $\sqrt{80} = 8.944$ , $\sqrt{40\times3.4142} = 11.686$ .
Pattern. In the lowest mode all three move in step, the middle mass most. In the middle mode the middle mass stays still while the outer two oppose. In the highest mode neighbours always oppose. The number of sign changes in the pattern is $0, 1, 2$ : higher frequency, more nodes. This is the picture of the harmonics of a string, with the masses as sample points.
The general formula. For $N$ masses the frequencies are $\omega_j = 2\sqrt{k/m}\,\sin\frac{j\pi}{2(N+1)}$ , $j = 1, \dots, N$ , with mode shapes $\big(\sin\frac{ji\pi}{N+1}\big)_{i=1}^N$ . For $N = 3$ these are $2\sqrt{40}\sin\frac{\pi}{8} = 4.841$ , $\ldots$ as above. Problem S2 proves this by one trigonometric identity. As $N \to \infty$ the highest frequency approaches $2\sqrt{k/m}$ and the lowest tends to $0$ like $\frac{1}{N+1}$ : the continuous string.

Answer. Frequencies $4.84$ , $8.94$ , $11.69$ rad/s, with shapes $(1, \sqrt2, 1)$ , $(1, 0, -1)$ , $(1, -\sqrt2, 1)$ .

Python labNormal-mode frequencies of a chain from the eigenvalues

import numpy as np

def chain_matrix(n, k=20.0, m=0.5):
    # n equal masses m joined by equal springs k, with a spring to a wall at each end
    M = np.zeros((n, n))
    for a in range(n):
        M[a, a] = 2 * k / m
        if a > 0:
            M[a, a - 1] = -k / m
        if a < n - 1:
            M[a, a + 1] = -k / m
    return M

k, m = 20.0, 0.5
print(" n    lowest w    highest w   (n+1)*lowest   largest gap to the formula")
for n in [1, 2, 3, 5, 10, 50, 200]:
    w = np.sqrt(np.linalg.eigvalsh(chain_matrix(n, k, m)))
    formula = 2 * np.sqrt(k / m) * np.sin(np.arange(1, n + 1) * np.pi / (2 * (n + 1)))
    print(f"{n:4d}  {w[0]:9.4f}  {w[-1]:10.4f}  {(n + 1) * w[0]:11.4f}   {np.max(np.abs(w - formula)):.1e}")
print("2*sqrt(k/m) =", 2 * np.sqrt(k / m), "  pi*sqrt(k/m) =", np.pi * np.sqrt(k / m))

Try this. Run it. The numbers in the last column are all below $10^{-13}$ : the eigenvalues of the matrix agree with the formula $2\sqrt{k/m}\sin\frac{j\pi}{2(N+1)}$ for every $N$ . The highest frequency creeps up to $2\sqrt{k/m} = 12.649$ and cannot pass it: a chain of masses cannot vibrate faster than that, whatever its length (the cut-off frequency of a lattice). The lowest frequency times $N + 1$ tends to $\pi\sqrt{k/m} = 19.869$ , which says the lowest frequency is about $19.87/(N+1)$ : a longer chain is softer. Change k or m and check the formula still holds, then change the line M[a, a] = 2 * k / m into M[a, a] = 3 * k / m (an extra spring from every mass to the ground) and watch the lowest frequency stop going to zero.

RLC circuits

A series circuit is a loop containing a resistor, an inductor, a capacitor and a voltage source. Kirchhoff’s voltage law says the source voltage equals the sum of the voltage drops around the loop. A resistor of $R$ ohms drops $RI$ when the current is $I$ . An inductor of $L$ henries drops $LI'$ : a changing current induces an opposing voltage. A capacitor of $C$ farads holding charge $q$ drops $q/C$ . The current is the rate of flow of charge, $I = q'$ .

Definition B

RLC circuit

The RLC circuit in series has the charge equation

L\,q'' + R\,q' + \frac1C\,q = E(t), \qquad I = q' .

It is the oscillator equation $mx'' + cx' + kx = F(t)$ with the dictionary

circuit	spring
charge $q$	displacement $x$
current $I = q'$	velocity $x'$
inductance $L$	mass $m$
resistance $R$	damping constant $c$
$1/C$	spring constant $k$
source voltage $E(t)$	force $F(t)$

The inductor stores energy $\frac12LI^2$ (kinetic energy), the capacitor $\frac{q^2}{2C}$ (potential energy), and the resistor turns energy into heat (friction).

In one lineA loop with a resistor $R$ , an inductor $L$ , a capacitor $C$ and a source $E(t)$ . The charge $q$ on the capacitor satisfies $Lq'' + Rq' + \dfrac qC = E(t)$ .

Why it existsIt is exactly the equation of a mass on a spring with friction and a driving force, so everything proved about oscillators in Unit 7 applies: damping, resonance, stability.

Write $\alpha = \frac{R}{2L}$ (the decay rate) and $\omega_0 = \frac{1}{\sqrt{LC}}$ (the natural frequency without damping). The characteristic polynomial $Lr^2 + Rr + \frac1C$ has roots $r = -\alpha \pm \sqrt{\alpha^2 - \omega_0^2}$ , so by 7.2 the free motion is underdamped (decaying oscillation $e^{-\alpha t}\cos\omega_dt$ with $\omega_d = \sqrt{\omega_0^2 - \alpha^2}$ ) when $R < 2\sqrt{L/C}$ , critically damped when $R = 2\sqrt{L/C}$ , and overdamped (two decaying exponentials) when $R > 2\sqrt{L/C}$ . In the language of 7.6: with $x = (q, I)$ the system is $x' = \begin{bmatrix}0 & 1\\-\frac{1}{LC} & -\frac RL\end{bmatrix}x + \begin{bmatrix}0\\E/L\end{bmatrix}$ with trace $-\frac RL < 0$ and determinant $\frac{1}{LC} > 0$ , so the origin is always asymptotically stable: a stable spiral when underdamped and a stable node when overdamped.

Fact B standard; Lyapunov function from 7.6

A passive circuit loses energy

Along a solution of $Lq'' + Rq' + q/C = E(t)$ with $I = q'$ , $\dot W = IE(t) - RI^2$ . Without a source, $W$ is a strict Lyapunov function away from $I = 0$ , and the circuit runs down.

In one lineWith no source, the stored energy $W = \frac12LI^2 + \frac{q^2}{2C}$ satisfies $\dot W = -RI^2 \le 0$ . With a source, $\dot W = IE - RI^2$ .

Why it is trueThe power the source delivers is $IE$ and the power dissipated in the resistor is $RI^2$ ; the rest is stored. As a function of the state it is a Lyapunov function.

Proofstandard

The idea: Differentiate W, replace LI' by the circuit equation.

\dot W = I\big(LI' + \frac qC\big)

Chain rule:

\dot W = LII' + \frac{qq'}{C}

, and

q' = I

\dot W = I(E - RI) = IE - RI^2

The circuit equation says

LI' + RI + \frac qC = E

, so

LI' + \frac qC = E - RI

\blacksquare

The forced circuit is what makes a radio work. Let the source be $E(t) = E_0\cos\omega t$ (a pure tone). Then the equation is linear with constant coefficients and a sinusoidal right-hand side, and 7.3 solves it by undetermined coefficients. The result is best written with complex numbers.

Theorem T 7.3 (undetermined coefficients, as in MATH 376 homework 9 #2); the complex impedance form was not in the course

Steady-state response of the driven RLC circuit

Let $L, R, C > 0$ and $\omega > 0$ . The current in $L q'' + Rq' + q/C = E_0\cos\omega t$ has the form $I(t) = I_{tr}(t) + I_{ss}(t)$ , where $I_{tr}(t) \to 0$ exponentially and $I_{ss}(t) = \operatorname{Re}\big(\frac{E_0}{Z}e^{i\omega t}\big) = \frac{E_0}{\lvert Z\rvert}\cos(\omega t - \varphi)$ with $\tan\varphi = \frac{\omega L - 1/(\omega C)}{R}$ and $\lvert Z\rvert = \sqrt{R^2 + \big(\omega L - \frac{1}{\omega C}\big)^2}$ .

In one lineFor $E(t) = E_0\cos\omega t$ , every solution is a decaying transient plus the steady state $I_{ss}(t) = \dfrac{E_0}{\lvert Z\rvert}\cos(\omega t - \varphi)$ , where $Z = R + i\left(\omega L - \dfrac{1}{\omega C}\right)$ is the impedance.

Why it is trueTry a complex exponential $\hat Ie^{i\omega t}$ in the equation for the current. Each derivative multiplies by $i\omega$ , so the equation becomes one complex equation $Z\hat I = E_0$ . The transient dies out because the characteristic roots have negative real part.

Proof7.3; standard

The idea: Differentiate the equation to get one for the current, try a complex exponential, solve for its amplitude, and use that the homogeneous solutions decay.

LI'' + RI' + \frac IC = E'(t) = -E_0\omega\sin\omega t = \operatorname{Re}\big(i\omega E_0e^{i\omega t}\big)

Differentiate

Lq'' + Rq' + q/C = E

once and use

I = q'

\left(L(i\omega)^2 + R(i\omega) + \frac1C\right)\hat I = i\omega E_0

Solve the complex equation

LI'' + RI' + I/C = i\omega E_0e^{i\omega t}

and take the real part (the coefficients are real, so the real part of a solution solves the real equation). Try

I = \hat Ie^{i\omega t}

; each derivative multiplies by

i\omega

(7.2). This is possible because

i\omega

is not a characteristic root (the roots have negative real part).

\left(R + i\omega L - \frac{i}{\omega C}\right)\hat I = E_0

, that is

Z\hat I = E_0

with

Z = R + i\big(\omega L - \frac1{\omega C}\big)

Divide both sides by

i\omega

\frac{L(i\omega)^2}{i\omega} = i\omega L

and

\frac{1}{C\,i\omega} = -\frac{i}{\omega C}

I_{ss}(t) = \frac{E_0}{\lvert Z\rvert}\cos(\omega t - \varphi)

is a solution.

Polar form of

Z = \lvert Z\rvert e^{i\varphi}

(0.6):

\hat I = \frac{E_0}{\lvert Z\rvert}e^{-i\varphi}

, and the real part of

\hat Ie^{i\omega t}

\frac{E_0}{\lvert Z\rvert}\cos(\omega t - \varphi)

The pair

(q_{ss}, I_{ss})

solves the circuit equation.

We solved the differentiated equation, so we still owe a charge. Define

q_{ss} = C\big(E - LI_{ss}' - RI_{ss}\big)

. Then

q_{ss}

satisfies the circuit equation by construction, and

q_{ss}' = C\big(E' - LI_{ss}'' - RI_{ss}'\big) = I_{ss}

because

I_{ss}

solves the differentiated equation. So

I_{ss}

really is the current of a solution.

The difference tends to

0

: that is the transient.

Any other solution differs from it by a solution of the homogeneous equation (7.2), which is a combination of

e^{r_1t}

e^{r_2t}

(or

e^{-\alpha t}\cos\omega_dt

etc.) with

\operatorname{Re}r_j = -\alpha

or smaller, negative.

\blacksquare

The current amplitude $\frac{E_0}{\lvert Z\rvert}$ is largest when the reactance $\omega L - \frac{1}{\omega C}$ vanishes, that is at $\omega = \omega_0 = \frac{1}{\sqrt{LC}}$ . There $\lvert Z\rvert = R$ and the amplitude is $\frac{E_0}{R}$ : resonance. How sharp the peak is depends on $R$ .

Definition B

Quality factor

The quality factor of the series RLC circuit is $Q = \frac{\omega_0L}{R}$ . The current amplitude falls to $\frac{1}{\sqrt2}$ of its peak (half the power) at the two frequencies $\omega$ with $\lvert\omega L - \frac1{\omega C}\rvert = R$ , and the distance between them is exactly $\Delta\omega = \frac RL$ . Hence $\frac{\Delta\omega}{\omega_0} = \frac1Q$ . A free oscillation has its amplitude multiplied by $e^{-1}$ in the time $\frac1\alpha = \frac{2L}{R}$ , during which the circuit completes about $\frac Q\pi$ oscillations.

In one line $Q = \dfrac{\omega_0L}{R} = \dfrac1R\sqrt{\dfrac LC}$ for a series RLC circuit: the ratio of the resonant frequency to the width of the resonance peak.

Why it existsIt measures how many oscillations the circuit makes before the energy dies out, and how selective a tuned circuit is.

Worked example MATH 376 homework 9 #2; 7.3

Homework 9 #2 of MATH 376 as a circuit: q″ + 5q′ + 4q = 10 sin 2t. Identify L, R, C and E, classify the damping, find the steady state, and explain the amplitude with the impedance.

Solution.

Read the circuit. $L = 1$ H, $R = 5\ \Omega$ , $\frac1C = 4$ so $C = 0.25$ F, and $E(t) = 10\sin2t$ V. Compare with $Lq'' + Rq' + q/C = E$ . (These numbers are far from those of a pocket radio, and the mathematics is the same.)
Classify the damping. $2\sqrt{L/C} = 2\sqrt{4} = 4 < R = 5$ : overdamped. The roots of $r^2 + 5r + 4 = (r+1)(r+4)$ are $-1$ and $-4$ . Two real negative roots: stable node in the phase plane.
The steady state, from 7.3. $q_p = -\cos2t$ , so the whole solution is $q = C_1e^{-t} + C_2e^{-4t} - \cos2t$ . The two exponentials are the transient and die out.
The steady-state current. $I_{ss} = q_p' = 2\sin2t$ , amplitude $2$ A. Differentiate the charge.
Check with the impedance. $\omega = 2$ : $\omega L = 2$ and $\frac{1}{\omega C} = \frac42 = 2$ . The reactance is $2 - 2 = 0$ , so $\lvert Z\rvert = R = 5$ and the amplitude is $\frac{E_0}{R} = \frac{10}{5} = 2$ . The circuit is driven exactly at its natural frequency $\omega_0 = \frac{1}{\sqrt{LC}} = 2$ : resonance, current in phase with the source ( $\varphi = 0$ ). The two methods agree.

Answer. $q(t) = C_1e^{-t} + C_2e^{-4t} - \cos2t$ ; steady-state current $2\sin2t$ , at resonance with amplitude $E_0/R$ .

Worked example the amplitude formula above

An AM radio tuner: L = 200 μH, R = 10 Ω, tuned to 1000 kHz. Find C, the quality factor, the bandwidth, and how well a station 10 kHz away is rejected.

Solution.

Find $C$ . $\omega_0 = 2\pi\times10^6$ rad/s, and $\omega_0^2 = \frac{1}{LC}$ gives $C = \frac{1}{\omega_0^2L} = \frac{1}{(2\pi\times10^6)^2\times2\times10^{-4}} = 126.65$ pF. Solve $\omega_0 = 1/\sqrt{LC}$ for $C$ .
Quality factor. $Q = \frac{\omega_0L}{R} = \frac{2\pi\times10^6\times2\times10^{-4}}{10} = 125.7$ . The definition.
Bandwidth. $\Delta f = \frac{f_0}{Q} = \frac{R}{2\pi L} = 7.96$ kHz. The half-power frequencies are $996.03$ kHz and $1003.99$ kHz. $\Delta\omega = R/L$ exactly; divide by $2\pi$ to get hertz.
The decay time. $\alpha = \frac{R}{2L} = 25\,000\ \text{s}^{-1}$ , so a free oscillation dies by the factor $e^{-1}$ in $\frac1\alpha = 40\ \mu$ s, which is $40$ cycles of $1$ MHz. Underdamped, since $Q \gg \frac12$ ; $\omega_d = \sqrt{\omega_0^2 - \alpha^2}$ is within $10^{-5}$ of $\omega_0$ .
A station 10 kHz away. At $f = 1.01$ MHz: $\omega L = 1269.2\ \Omega$ and $\frac{1}{\omega C} = 1244.2\ \Omega$ , so the reactance is $25.0\ \Omega$ and $\lvert Z\rvert = \sqrt{10^2 + 25.0^2} = 26.9\ \Omega$ . The current is $\frac{10}{26.9} = 0.37$ of what the same source would produce at resonance (about $-8.6$ dB). The formula for $\lvert Z\rvert$ at a detuned frequency.

Answer. $C = 126.7$ pF, $Q = 125.7$ , bandwidth $7.96$ kHz, a neighbouring channel is cut to $37\%$ by this one stage. A radio uses several such stages.

Your answerFor the tuner above (

L = 200\ \mu

R = 10\ \Omega

f_0 = 1

MHz), what is the quality factor

Q = \omega_0L/R

? (Two decimals.)

Your answerIn the tuner, at

1.01

MHz the impedance is

26.9\ \Omega

while at resonance it is

10\ \Omega

. What fraction of the resonant current flows at

1.01

MHz? (Two decimals.)

Quick checkAn RLC circuit has

L = 0.1

C = 100\ \mu

F and

R = 20\ \Omega

. Is it underdamped or overdamped?

The SIR epidemic

Divide a population into three groups: Susceptible (can catch the disease), Infected (and infectious), Recovered (immune, or removed). Let $s$ , $i$ and $r$ be their fractions of the population, so $s + i + r = 1$ .

Definition B

SIR model

The SIR model of Kermack and McKendrick is

s' = -\beta\,s\,i, \qquad i' = \beta\,s\,i - \gamma\,i, \qquad r' = \gamma\,i,

with constants $\beta > 0$ (the infection rate: new infections per infected person per unit time in a fully susceptible population) and $\gamma > 0$ (the recovery rate: $1/\gamma$ is the mean time a person stays infectious). The sum $s + i + r$ is constant, because the three right-hand sides add to $0$ , so the model is really the planar system for $(s, i)$ .

In one lineThe system $s' = -\beta si$ , $i' = \beta si - \gamma i$ , $r' = \gamma i$ for the fractions of susceptible, infected and recovered people.

Why it existsIt is the simplest model of an epidemic that keeps track of who can still be infected. Infection needs a meeting between an infected and a susceptible person (the product $si$ ), recovery takes a fixed fraction of the infected per unit time.

Definition B

Basic reproduction number

The basic reproduction number is $R_0 = \dfrac\beta\gamma$ .

In one line $R_0 = \beta/\gamma$ : the average number of people infected by one infected person during the whole infectious period, when everyone else is susceptible.

Why it existsOne infected person infects $\beta$ people per unit time, for a time $1/\gamma$ on average. The product is the number of 'children' of an infection, and whether it is above or below $1$ decides whether the infection grows.

Definition B

Disease-free equilibrium

A disease-free equilibrium of the SIR model is a point $(s^*, i^*) = (s^*, 0)$ with $0 \le s^* \le 1$ . Each is an equilibrium, since $s' = -\beta s\cdot0 = 0$ and $i' = (\beta s - \gamma)\cdot0 = 0$ . (There are no others: $i' = 0$ needs $i = 0$ or $s = 1/R_0$ , and $s' = 0$ then forces $i = 0$ .)

In one lineA state with no infected people, $i = 0$ . In the SIR model every $(s^*, 0)$ is an equilibrium, so these form a whole line.

Why it existsIt is the state the epidemic starts from and ends in. The question is whether it is stable: does a few infected people grow into an epidemic?

Theorem B Kermack and McKendrick 1927; linearization (7.6)

The epidemic threshold

For the SIR system in $(s, i)$ the Jacobian at $(s^*, 0)$ is $\begin{bmatrix}0 & -\beta s^*\\0 & \beta s^* - \gamma\end{bmatrix}$ with eigenvalues $0$ and $\lambda = \beta s^* - \gamma = \gamma(R_0s^* - 1)$ . If $s^* > 1/R_0$ then $\lambda > 0$ and a small infection grows like $e^{\lambda t}$ ; if $s^* < 1/R_0$ it decays. In a fully susceptible population ( $s^* = 1$ ) there is an epidemic exactly when $R_0 > 1$ .

In one lineAt the disease-free equilibrium $(s^*, 0)$ the Jacobian of the $(s, i)$ system has eigenvalues $0$ and $\gamma(R_0s^* - 1)$ . A few infected people grow exponentially if $R_0s^* > 1$ and die out if $R_0s^* < 1$ .

Why it is trueThe eigenvalue $0$ is the line of equilibria itself. The other eigenvalue is the growth rate of the infection: with $s \approx s^*$ the equation $i' = (\beta s - \gamma)i$ is a linear equation with growth rate $\beta s^* - \gamma$ .

Proof7.6

The idea: Differentiate the two right-hand sides to get the Jacobian. At i = 0 it is triangular, so the eigenvalues are the diagonal entries.

DF = \begin{bmatrix}-\beta i & -\beta s\\\beta i & \beta s - \gamma\end{bmatrix}

Rows are the gradients of

-\beta si

and

\beta si - \gamma i

(5.3):

\frac{\partial}{\partial s}(-\beta si) = -\beta i

\frac{\partial}{\partial i}(-\beta si) = -\beta s

, and so on.

DF(s^*, 0) = \begin{bmatrix}0 & -\beta s^*\\0 & \beta s^* - \gamma\end{bmatrix}

Substitute

i = 0

s = s^*

. The first column becomes

0

The eigenvalues are

0

and

\beta s^* - \gamma = \gamma(R_0s^* - 1)

, because

\beta = \gamma R_0

The matrix is upper triangular: the eigenvalues are the diagonal entries (2.8). The eigenvector for

0

(1, 0)

, pointing along the line of equilibria.

If the second eigenvalue is positive the infection grows; if negative it decays.

Near the equilibrium

i' = (\beta s - \gamma)i

and

s \approx s^*

varies slowly while

i

is small (

s' = -\beta si

is of the size of

i

). So

i

grows like

e^{(\beta s^* - \gamma)t}

if the exponent is positive. The equilibrium is not hyperbolic (the eigenvalue

0

is forced by the line of equilibria), so this conclusion is about the transverse direction.

\blacksquare

For $R_0 = 3$ and $\gamma = 0.2$ per day (infectious for five days on average, $\beta = 0.6$ per day) the early growth rate is $\gamma(R_0 - 1) = 0.4$ per day, and the number of cases doubles every $\frac{\ln2}{0.4} = 1.73$ days. A fraction $p$ of immunity at the start gives $s^* = 1 - p$ , so the infection dies out as soon as $R_0(1 - p) < 1$ , that is $p > 1 - \frac{1}{R_0}$ : the herd immunity threshold. For $R_0 = 3$ it is $\frac23$ .

Worked example Kermack and McKendrick 1927; new numbers

An infection has β = 0.35 per day and γ = 0.1 per day. Find R0 and the herd immunity threshold. Does an epidemic start if 50% of the population is immune, and if 80% is? Give the growth rate or the decay rate in each case.

Solution.

Compute $R_0$ . $R_0 = \frac{\beta}{\gamma} = \frac{0.35}{0.1} = 3.5$ . Definition of the basic reproduction number. One infected person infects $0.35$ people per day for an average of $10$ days.
Herd immunity threshold. $1 - \frac{1}{R_0} = 1 - \frac{1}{3.5} = 0.714$ , that is $71.4\%$ . The infection dies out as soon as $R_0(1 - p) < 1$ .
Fifty percent immune. Then $s^* = 0.5$ and $R_0s^* = 1.75 > 1$ . The growth eigenvalue is $\gamma(R_0s^* - 1) = 0.1\times0.75 = 0.075$ per day, so the number of cases doubles every $\frac{\ln2}{0.075} = 9.2$ days. The second eigenvalue of the Jacobian at $(s^*, 0)$ is positive: a small infection grows.
Eighty percent immune. Then $s^* = 0.2$ and $R_0s^* = 0.7 < 1$ . The eigenvalue is $0.1\times(0.7 - 1) = -0.03$ per day, so a small infection shrinks by the factor $e^{-0.03t}$ , halving every $\frac{\ln2}{0.03} = 23$ days. The eigenvalue is negative, and $80\%$ is above the threshold of $71.4\%$ , as it must be.

Answer. $R_0 = 3.5$ ; threshold $71.4\%$ ; at $50\%$ immune the infection grows at $0.075$ per day (doubling in $9.2$ days); at $80\%$ immune it decays at $0.03$ per day.

The linearization says how the epidemic starts. The second theorem says how it develops and ends, by a conserved quantity.

Theorem B Kermack and McKendrick 1927; Murray 10

Peak and final size of an epidemic

Let $s_0, i_0 > 0$ be the starting values. (a) $\dot K = 0$ for $K = i + s - \frac1{R_0}\ln s$ . (b) If $R_0s_0 > 1$ the fraction infected rises to a maximum at the moment $s = \frac1{R_0}$ and equals $i_{\max} = i_0 + s_0 - \frac{1 + \ln(R_0s_0)}{R_0}$ . (c) As $t \to \infty$ , $i \to 0$ and $s \to s_\infty > 0$ , the root in $(0, \frac1{R_0})$ of $\ln\frac{s_0}{s} = R_0(s_0 + i_0 - s)$ . In particular the epidemic ends before everyone has been infected.

In one lineAlong every solution of the SIR system, $i + s - \dfrac1{R_0}\ln s$ is constant. Hence the infection peaks when $s = 1/R_0$ at $i_{\max} = i_0 + s_0 - \dfrac{1 + \ln(R_0s_0)}{R_0}$ , and the epidemic ends with $i \to 0$ and $s \to s_\infty$ where $\ln\dfrac{s_0}{s_\infty} = R_0(s_0 + i_0 - s_\infty)$ .

Why it is trueDividing the two equations removes time: $\frac{di}{ds} = -1 + \frac{1}{R_0s}$ , which integrates at once. The conserved quantity ties the start of the epidemic to its end and its peak.

ProofKermack and McKendrick 1927

The idea: Differentiate K along solutions and watch the terms cancel. Evaluate K at the peak and at the end.

\dot K = -\gamma i + \gamma i = 0

Chain rule:

\dot K = i' + s' - \frac{s'}{R_0s}

. Insert

i' + s' = -\gamma i

(add the first two equations) and

\frac{s'}{R_0s} = \frac{-\beta si}{R_0s} = -\gamma i

, because

\frac\beta{R_0} = \gamma

The maximum of

i

is where

s = \frac1{R_0}

i' = i(\beta s - \gamma)

, which is zero when

s = \gamma/\beta = 1/R_0

(and

i > 0

). Since

s

decreases,

i

increases while

s > 1/R_0

and decreases afterwards.

i_{\max} = i_0 + s_0 - \frac{1}{R_0} - \frac{1}{R_0}\ln(R_0s_0) = i_0 + s_0 - \frac{1 + \ln(R_0s_0)}{R_0}

Use

K(s = 1/R_0, i_{\max}) = K(s_0, i_0)

i_{\max} + \frac1{R_0} - \frac{1}{R_0}\ln\frac1{R_0} = i_0 + s_0 - \frac1{R_0}\ln s_0

. Solve for

i_{\max}

At the end,

\ln\frac{s_0}{s_\infty} = R_0(s_0 + i_0 - s_\infty)

s

is decreasing and bounded below by

0

, so it has a limit

s_\infty

. Since

r' = \gamma i

and

r \le 1

, the integral

\int_0^\infty i\,dt = (r(\infty) - r(0))/\gamma

is finite. A non-negative function with finite integral and a bounded derivative (

\lvert i' \rvert \le (\beta + \gamma)

) tends to

0

, so

i \to 0

. Then

K(s_\infty, 0) = K(s_0, i_0)

gives

s_\infty - \frac{1}{R_0}\ln s_\infty = i_0 + s_0 - \frac1{R_0}\ln s_0

Hence

0 < s_\infty < \frac1{R_0}

: some people are never infected, and the epidemic overshoots the herd immunity level.

s_\infty

were

0

, then

-\frac1{R_0}\ln s \to +\infty

while

i, s \ge 0

, so

K \to +\infty

, but

K

is constant. And

s_\infty \le \frac{1}{R_0}

: if

s_\infty > \frac1{R_0}

, then

i' = i(\beta s - \gamma) > 0

for large

t

, so

i

would increase, contradicting

i \to 0

with

i > 0

. Equality

s_\infty = \frac1{R_0}

is also impossible: the function

s - \frac{1}{R_0}\ln s

has its minimum value

\frac{1 + \ln R_0}{R_0}

exactly at

s = \frac1{R_0}

, while

K(s_0, i_0)

is larger than that minimum because

i_0 > 0

(and

K(s_\infty, 0)

must equal

K(s_0, i_0)

\blacksquare

Proof workshopRebuild: the SIR quantity K is conservedKermack and McKendrick 1927

Each line of the proof is shown in order. Choose the reason that justifies it, then press Check.

$\dot K = i' + s' - \dfrac{s'}{R_0\,s}$ .
$i' + s' = -\gamma i$ .
$\dfrac{s'}{R_0s} = \dfrac{-\beta si}{R_0s} = -\dfrac{\beta}{R_0}\,i = -\gamma i$ .
$\dot K = -\gamma i - (-\gamma i) = 0$ .

Worked example Kermack and McKendrick 1927; new numbers

The SIR epidemic with R0 = 3: peak, final size and timing, then a table for other R0.

Solution. Take $\gamma = 0.2$ per day, $\beta = 0.6$ per day, so $R_0 = 3$ . Start with $i_0 = 0.001$ (one in a thousand infected) and $s_0 = 0.999$ .

Does an epidemic start? $R_0s_0 = 3\times0.999 = 2.997 > 1$ , eigenvalue $\gamma(R_0s_0 - 1) = 0.2\times1.997 = 0.399$ per day. Positive, so the infection grows at about $0.4$ per day, doubling in $1.7$ days.
The peak. $i_{\max} = 0.001 + 0.999 - \frac{1 + \ln(2.997)}{3} = 1 - \frac{1 + 1.0976}{3} = 0.3008$ . It occurs when $s = \frac13$ . Peak formula; at the peak a third of the population is still susceptible.
The final size. Solve $\ln\frac{0.999}{s_\infty} = 3(1 - s_\infty)$ . By iteration of $s \mapsto 0.999\,e^{-3(1 - s)}$ starting from $0.5$ : $0.5 \to 0.223 \to 0.097 \to 0.067 \to 0.061 \to \cdots$ , settling at $s_\infty = 0.0594$ . Fixed-point iteration for the root; it converges because the derivative $3s_\infty = 0.18$ is below $1$ .
Interpret. $94.1\%$ of the population is infected at some point, and only $5.9\%$ escapes. The herd immunity threshold is $\frac23$ but the epidemic overshoots: it ends at $s_\infty = 0.059$ , well below $\frac13$ . After the peak the infection still has many infectives, who keep infecting until $i$ has fallen to $0$ .
Timing. The numerical solution (lab below) puts the peak at $3.84/\gamma = 19.2$ days after the start, from $0.1\%$ infected. Time is measured in units of $1/\gamma = 5$ days.

$R_0$	herd immunity $1 - 1/R_0$	peak fraction infected	final fraction ever infected
$0.8$	none (no epidemic)	$0.001$ (the start)	$0.5\%$
$1.5$	$33\%$	$6.4\%$	$58.4\%$
$2$	$50\%$	$15.4\%$	$79.7\%$
$3$	$67\%$	$30.1\%$	$94.1\%$
$15$	$93\%$	$75.3\%$	$\approx100\%$

Answer. Peak $30.1\%$ infected at $19$ days, final $94.1\%$ ever infected, for $R_0 = 3$ .

Your answerFor the SIR model with

R_0 = 3

started from a fully susceptible population (

s_0 \approx 1

i_0 \approx 0

), what is the largest fraction of the population infected at one time? (Three decimals.)

Your answerMeasles has

R_0

of about

15

. What fraction of the population must be immune to stop an epidemic? (A fraction.)

Quick checkAn SIR epidemic has

R_0 = 2

. A vaccination campaign immunises 40% of the population before the disease arrives. Will an epidemic start?

Python labSimulating the SIR epidemic and checking the formulas

import numpy as np

I0 = 0.001            # fraction infected at the start

def sir(R0, gamma=1.0, T=60.0, h=0.01):
    # fractions s (susceptible) and i (infected); recovered r = 1 - s - i. Time is measured in units of 1/gamma.
    beta = R0 * gamma
    f = lambda s, i: (-beta * s * i, beta * s * i - gamma * i)
    s, i, t = 1.0 - I0, I0, 0.0
    peak, t_peak = i, 0.0
    for _ in range(int(T / h)):
        k1 = f(s, i)
        k2 = f(s + h/2*k1[0], i + h/2*k1[1])
        k3 = f(s + h/2*k2[0], i + h/2*k2[1])
        k4 = f(s + h*k3[0], i + h*k3[1])
        s += h/6*(k1[0] + 2*k2[0] + 2*k3[0] + k4[0])
        i += h/6*(k1[1] + 2*k2[1] + 2*k3[1] + k4[1])
        t += h
        if i > peak:
            peak, t_peak = i, t
    return s, peak, t_peak

def final_size(R0):
    # solve ln(s0 / s) = R0 (s0 + i0 - s) = R0 (1 - s) by the iteration s = s0 exp(-R0 (1 - s))
    s = 0.5
    for _ in range(500):
        s = (1 - I0) * np.exp(-R0 * (1 - s))
    return s

print(" R0   final s (sim)  final s (theory)  ever infected   peak i (sim)  peak i (theory)  peak time")
for R0 in [0.8, 1.5, 2.0, 3.0, 15.0]:
    s, peak, tp = sir(R0)
    th_peak = 1 - (1 + np.log(R0 * (1 - I0))) / R0 if R0 > 1 else I0
    print(f"{R0:4.1f}   {s:11.6f}   {final_size(R0):13.6f}   {1 - s:11.1%}   {peak:10.4f}   {th_peak:13.4f}   {tp:8.2f}")

Try this. Run it. For each $R_0$ the simulated final susceptible fraction agrees with the root of $\ln(s_0/s) = R_0(1 - s)$ to 6 digits, and the simulated peak agrees with $1 - (1 + \ln(R_0s_0))/R_0$ to 3 or 4 digits (the last digit at $R_0 = 15$ differs because of the step size). For $R_0 = 0.8$ there is no epidemic: the infected fraction only falls, and the total ever infected is $0.5\%$ , the starting $0.1\%$ plus a few more. For $R_0 = 15$ only about $3 \times 10^{-7}$ of the population escapes. Now change I0 = 0.001 to I0 = 0.0001: the peak and the final size hardly change ( $R_0 = 3$ still gives a peak of $30.0\%$ and $94.0\%$ ever infected), but the peak comes later, at $4.99$ instead of $3.84$ in units of $1/\gamma$ , because the infection takes longer to grow from a smaller start. Restore I0 and change T=60.0 to T=10.0: for $R_0 = 1.5$ the epidemic is still in progress at time $10$ (the susceptible fraction is $0.758$ , not $0.416$ ).

A two-level quantum system

In quantum mechanics the state of a system is a vector $\psi$ of unit length in a complex vector space, and the energy is an operator. For a system with only two relevant states (the two spin directions of an electron or proton, the ground and excited state of an atom driven by a laser, the two lowest states of a superconducting circuit) the state space is $\mathbb C^2$ and the operator is a $2 \times 2$ matrix.

Definition B

Hamiltonian

The Hamiltonian is a Hermitian matrix $H$ (3.1), $H^* = H$ , so its eigenvalues $E_1, E_2, \dots$ are real: the allowed energies. The state $\psi(t)$ obeys the Schrödinger equation

i\hbar\,\frac{d\psi}{dt} = H\psi, \qquad \hbar = 1.0546\times10^{-34}\ \text{J s} .

The probability of finding the system in the state $u$ (a unit vector) is $\lvert u^*\psi\rvert^2$ .

In one lineThe Hermitian matrix $H$ (or operator) of a quantum system whose eigenvalues are the possible energies. The state obeys the Schrödinger equation $i\hbar\,\psi'(t) = H\psi(t)$ .

Why it existsIt packages the physics of the system, the forces and fields, into one matrix. Everything about the motion follows from it by the matrix exponential.

Definition B

Two-level system

A two-level system has a state $\psi \in \mathbb C^2$ with $\lvert\psi_1\rvert^2 + \lvert\psi_2\rvert^2 = 1$ ( $\lvert\psi_1\rvert^2$ is the probability of level 1) and a $2 \times 2$ Hermitian matrix $H$ . The Pauli matrices $\sigma_x = \begin{bmatrix}0 & 1\\1 & 0\end{bmatrix}$ and $\sigma_z = \begin{bmatrix}1 & 0\\0 & -1\end{bmatrix}$ appear in nearly every Hamiltonian.

In one lineA quantum system whose state space is $\mathbb C^2$ : a vector $\psi = (\psi_1, \psi_2)$ with $\lvert\psi_1\rvert^2 + \lvert\psi_2\rvert^2 = 1$ , and a $2 \times 2$ Hermitian Hamiltonian.

Why it existsIt is the smallest quantum system, and the model of a spin- $\frac12$ particle in a magnetic field, of a qubit, and of an atom driven near a single resonance.

The Schrödinger equation is $\psi' = A\psi$ with the constant matrix $A = -\frac{i}{\hbar}H$ (a complex matrix; everything in 7.4 works for complex matrices, as MATH 376’s Remark 8.10 says). So $\psi(t) = e^{-iHt/\hbar}\psi(0)$ .

Theorem B standard; 7.4 and 3.3

Unitary evolution

Let $H$ be a Hermitian $n \times n$ matrix. Then $\psi(t) = e^{-iHt/\hbar}\psi(0)$ is the unique solution of $i\hbar\psi' = H\psi$ , the matrix $U(t) = e^{-iHt/\hbar}$ is unitary (3.3) for every real $t$ , and $\psi$ keeps its length. If $H = V\operatorname{diag}(E_j)V^*$ with $V$ unitary, then $U(t) = V\operatorname{diag}\big(e^{-iE_jt/\hbar}\big)V^*$ : each energy eigenstate is only multiplied by a phase.

In one lineThe solution of $i\hbar\psi' = H\psi$ is $\psi(t) = U(t)\psi(0)$ with $U(t) = e^{-iHt/\hbar}$ . For Hermitian $H$ the matrix $U(t)$ is unitary, so $\lVert\psi(t)\rVert = \lVert\psi(0)\rVert$ for all $t$ : total probability stays $1$ .

Why it is trueThe adjoint of the exponential series is the series of the adjoint, and $H^* = H$ turns $-iHt/\hbar$ into $+iHt/\hbar$ . So $U^* = e^{+iHt/\hbar} = U^{-1}$ .

Proofstandard

The idea: Solve with the matrix exponential, take the adjoint of the series, and use the law of exponents for the commuting matrices iHt/hbar and -iHt/hbar.

\psi(t) = e^{-iHt/\hbar}\psi(0)

The equation is

\psi' = A\psi

with

A = -\frac i\hbar H

constant, and 7.4 gives

\psi(t) = e^{tA}\psi(0)

(unique by 7.1).

U^* = e^{iHt/\hbar}

Conjugate transpose each term of the series:

\big(\frac{1}{k!}(-iHt/\hbar)^k\big)^* = \frac1{k!}(iH^*t/\hbar)^k

because

(cB)^* = \bar cB^*

and

(B^k)^* = (B^*)^k

. Then use

H^* = H

. The adjoint is continuous, so it passes through the sum.

U^*U = e^{iHt/\hbar}e^{-iHt/\hbar} = I

, so

U

is unitary.

iHt/\hbar

and

-iHt/\hbar

commute (one is minus the other), so

e^{A}e^{-A} = e^0 = I

(7.4).

\lVert\psi(t)\rVert = \lVert\psi(0)\rVert = 1

A unitary matrix preserves lengths:

\lVert U\psi\rVert^2 = \psi^*U^*U\psi = \psi^*\psi

. For the diagonal form:

H = VDV^*

gives

H^k = VD^kV^*

and the series sums to

V\operatorname{diag}(e^{-iE_jt/\hbar})V^*

\blacksquare

In the language of 7.6: the eigenvalues of $-\frac{i}{\hbar}H$ are $-\frac{iE_j}{\hbar}$ , purely imaginary, the signature of a centre. The state rotates and neither decays nor grows. Real systems lose phase coherence to their surroundings, which adds damping and is the main enemy of quantum computers.

Theorem B Rabi 1937; Feynman Lectures III, chapters 8 to 11

The Rabi formula

Let $N = \delta\sigma_z + \Omega\sigma_x = \begin{bmatrix}\delta & \Omega\\\Omega & -\delta\end{bmatrix}$ with real $\delta, \Omega$ and $\Omega_R = \sqrt{\delta^2 + \Omega^2} > 0$ , and $H = \frac\hbar2N$ . Then $U(t) = e^{-iHt/\hbar} = \cos\frac{\Omega_Rt}2\,I - i\sin\frac{\Omega_Rt}2\,\frac{N}{\Omega_R}$ , and for $\psi(0) = (1, 0)$ the probability of finding the system in level 2 at time $t$ is $P(t) = \frac{\Omega^2}{\Omega_R^2}\sin^2\frac{\Omega_Rt}{2}$ . It oscillates between $0$ and $P_{\max} = \frac{\Omega^2}{\delta^2 + \Omega^2}$ with angular frequency $\Omega_R$ .

In one lineFor $H = \frac\hbar2(\delta\sigma_z + \Omega\sigma_x)$ and $\Omega_R = \sqrt{\delta^2 + \Omega^2}$ : $e^{-iHt/\hbar} = \cos\frac{\Omega_Rt}{2}\,I - i\sin\frac{\Omega_Rt}{2}\cdot\frac{\delta\sigma_z + \Omega\sigma_x}{\Omega_R}$ . Starting in level 1, the probability of level 2 is $\dfrac{\Omega^2}{\Omega_R^2}\sin^2\dfrac{\Omega_Rt}{2}$ .

Why it is trueThe matrix $N = \delta\sigma_z + \Omega\sigma_x$ satisfies $N^2 = \Omega_R^2I$ , exactly the situation of the $2\times2$ matrix exponential of 7.4 (a square that is a multiple of $I$ ), so the series splits into cosine and sine.

Proof7.4

The idea: Square N by hand, then sum the exponential series in even and odd powers.

N^2 = \begin{bmatrix}\delta^2 + \Omega^2 & 0\\0 & \Omega^2 + \delta^2\end{bmatrix} = \Omega_R^2I

Direct multiplication: the diagonal entries are

\delta^2 + \Omega^2

and the off-diagonal entries are

\delta\Omega - \Omega\delta = 0

e^{-iHt/\hbar} = \cos s\,I - i\sin s\,\frac{N}{\Omega_R}

with

s = \frac{\Omega_Rt}{2}

Put

B = -i\frac t2N = -iHt/\hbar

and

s = \frac{\Omega_Rt}{2}

. Then

B^2 = -\frac{t^2}{4}\Omega_R^2I = -s^2I

, so

B^{2k} = (-1)^ks^{2k}I

and

B^{2k+1} = (-1)^ks^{2k}B

. The series

\sum\frac{B^n}{n!}

splits into even and odd powers, exactly as for the rotation matrix in 7.4:

e^B = \cos s\,I + \frac{\sin s}{s}B

. And

\frac Bs = -i\frac{(t/2)N}{\Omega_Rt/2} = -i\frac{N}{\Omega_R}

\psi(t) = \Big(\cos s - i\frac{\delta}{\Omega_R}\sin s,\ -i\frac{\Omega}{\Omega_R}\sin s\Big)

Apply to

(1, 0)^{\mathsf T}

: the first column of

N

(\delta, \Omega)^{\mathsf T}

P(t) = \lvert\psi_2\rvert^2 = \frac{\Omega^2}{\Omega_R^2}\sin^2s

The squared modulus of the second component. As a check,

\cos^2s + \frac{\delta^2}{\Omega_R^2}\sin^2s + \frac{\Omega^2}{\Omega_R^2}\sin^2s = 1

\blacksquare

Worked example Rabi 1937; new numbers

A spin-1/2 or qubit driven by a field: find the pulse that flips it, the effect of detuning, and the numbers for a proton in an MRI magnet.

Solution.

On resonance, $\delta = 0$ . Then $\Omega_R = \Omega$ and $P(t) = \sin^2\frac{\Omega t}{2}$ . The probability is $1$ when $\Omega t = \pi$ : a $\pi$ -pulse turns level 1 completely into level 2. At $\Omega t = \frac\pi2$ the probability is $\sin^2\frac\pi4 = \frac12$ : an equal superposition (a $\frac\pi2$ -pulse). Set $\frac{\Omega t}{2} = \frac\pi2$ and $\frac\pi4$ .
Numbers for a qubit. With a Rabi frequency $\Omega/2\pi = 10$ MHz, $\Omega = 6.283\times10^7$ rad/s, the $\pi$ -pulse lasts $t = \frac\pi\Omega = 50$ ns and the $\frac\pi2$ -pulse $25$ ns. Divide by $\Omega$ . Superconducting qubits are driven at rates of this order, which is why a gate takes tens of nanoseconds.
Off resonance. If the drive is detuned by $\delta = \Omega$ then $\Omega_R = \sqrt2\,\Omega$ and $P_{\max} = \frac{\Omega^2}{2\Omega^2} = \frac12$ : the flip is never complete. For $\delta = 3\Omega$ it is only $\frac{1}{10}$ . The Rabi formula: a drive must be tuned to the transition. This is also why the quantum circuit behaves as a narrow filter, like the RLC circuit.
The MRI proton. A proton in a magnetic field $B_0 = 1.5$ T has two energy levels split by $\hbar\omega_L$ with $\frac{\omega_L}{2\pi} = 42.577\ \text{MHz/T}\times1.5\ \text{T} = 63.87$ MHz (the Larmor frequency), so the Hamiltonian is $H = \frac{\hbar\omega_L}{2}\sigma_z$ and $e^{-iHt/\hbar} = \operatorname{diag}(e^{-i\omega_Lt/2}, e^{i\omega_Lt/2})$ : a pure rotation of phase. The diagonal case of the unitary evolution theorem.
The radio-frequency pulse. A rotating field of amplitude $B_1 = 10\ \mu$ T at the Larmor frequency drives the spin at the Rabi frequency $\frac{\Omega}{2\pi} = 42.577\ \text{MHz/T}\times10^{-5}\ \text{T} = 425.8$ Hz. A $90^\circ$ pulse ( $\Omega t = \frac\pi2$ ) lasts $\frac{1}{4\times425.8\ \text{Hz}} = 587\ \mu$ s. $t = \frac{\pi/2}{\Omega} = \frac{1}{4(\Omega/2\pi)}$ . These are the numbers that set the timing of MRI sequences.

Answer. $P(t) = \sin^2\frac{\Omega t}{2}$ on resonance; $50$ ns for a $\pi$ -pulse at $10$ MHz; $\frac{\Omega^2}{\delta^2 + \Omega^2}$ off resonance; $587\ \mu$ s for a $90^\circ$ pulse of $10\ \mu$ T on a proton.

Your answerA qubit is driven on resonance with Rabi frequency

\Omega/2\pi = 10

MHz. How long, in nanoseconds, is a

\pi

-pulse (

\Omega t = \pi

Quick checkFor

\psi(t) = e^{-iHt/\hbar}\psi(0)

with

H

Hermitian, which quantity does NOT change with time?

Where it is used

Used inside maths· Inside maths: a network of springs is a graph Laplacian, and energy is the Lyapunov function

Join any number of unit masses by unit springs along the edges of a graph. The stiffness matrix is the graph Laplacian $\mathcal L = D - A_G$ (degrees on the diagonal, minus the adjacency matrix), and the normal modes are its eigenvectors. For the three masses in a row with no walls, $\mathcal L = \begin{bmatrix}1 & -1 & 0\\-1 & 2 & -1\\0 & -1 & 1\end{bmatrix}$ has eigenvalues $0, 1, 3$ with modes $(1, 1, 1)$ (the whole chain translating: frequency $0$ ), $(1, 0, -1)$ and $(1, -2, 1)$ . Compare with the chain between walls above: the walls remove the zero mode and raise the frequencies. In data science, the lowest non-zero mode (the Fiedler vector, here $(1, 0, -1)$ ) splits a graph into two halves, which is spectral clustering.

Energy plays the role of Lyapunov function in all four systems of this chapter. For the springs the energy is conserved (a centre: stable, not asymptotically). For the RLC circuit $\dot W = -RI^2$ (a stable spiral or node). For the SIR model $K$ is conserved and decides the final size. For the quantum state $\lVert\psi\rVert^2$ is conserved. In each case the conserved or decreasing quantity is the shortest route to the qualitative behaviour, and the eigenvalues give the quantitative detail.

Used in the real world· Civil engineering: the modes of a two-storey building

Model a two-storey building as two floors of mass $m = 10^5$ kg (100 tonnes), each connected to the one below by columns of stiffness $k = 4\times10^7$ N/m (illustrative values). Then $m\,x'' = -Kx$ with $K = k\begin{bmatrix}2 & -1\\-1 & 1\end{bmatrix}$ : the lower floor is held by the ground floor and the floor above, the upper floor only by the floor below. With $\frac km = 400\ \text{s}^{-2}$ the eigenvalues of $\frac Km$ are $400\cdot\frac{3 \mp \sqrt5}{2}$ , so $\omega_1 = 12.36$ rad/s ( $1.97$ Hz) and $\omega_2 = 32.36$ rad/s ( $5.15$ Hz). In the first mode the upper floor moves $1.618$ times as far as the lower one, in the same direction. In the second the floors move in opposite directions. The ratio of the two eigenvalues is $\frac{3 + \sqrt5}{3 - \sqrt5} = 6.854$ . Ground motion of an earthquake contains frequencies from below $1$ Hz to $10$ Hz, so the building is shaken hardest when the ground has energy near $2$ Hz or $5$ Hz; the design task is to move these frequencies away from the dominant frequencies of the expected ground motion, or to add damping so that the $Q$ of each mode is small.

Used in the real world· Radio: covering the AM band with one coil

An AM receiver tunes with a fixed inductor and a variable capacitor. The band runs from about $540$ kHz to $1600$ kHz. Since $f_0 = \frac{1}{2\pi\sqrt{LC}}$ , the capacitance must change by the factor $\left(\frac{1600}{540}\right)^2 = 8.78$ . With $L = 240\ \mu$ H the capacitor must cover $C = \frac{1}{(2\pi f)^2L}$ , that is $362$ pF at $540$ kHz down to $41$ pF at $1600$ kHz. With the series resistance $R = 10\ \Omega$ of the coil the quality factor at $1000$ kHz is $Q = \frac{2\pi\times10^6\times2.4\times10^{-4}}{10} = 151$ and the bandwidth is $\frac{f_0}{Q} = 6.6$ kHz, narrower than the $10$ kHz spacing of AM stations in the Americas, so that a neighbouring station is attenuated.

Used in the real world· Public health: why measles needs about 95% coverage

Measles has a reproduction number $R_0$ that is often quoted as between $12$ and $18$ (estimates vary with the place and the time). The herd immunity threshold $1 - \frac{1}{R_0}$ is then between $91.7\%$ and $94.4\%$ : with $R_0 = 15$ it is $93.3\%$ . Because a vaccine does not protect everyone who receives it, the proportion that has to be vaccinated is higher than the proportion that has to be immune; this is why the public-health target for two-dose measles vaccination is about $95\%$ . The threshold theorem is also a statement about the early phase. Take $R_0 = 15$ and a community with $90\%$ immune: then $s^* = 0.1$ and the growth eigenvalue is $\gamma(R_0s^* - 1) = \gamma\times0.5 > 0$ , so an outbreak starts and grows, slowly. At $95\%$ immune the eigenvalue is $\gamma(15\times0.05 - 1) = -0.25\gamma < 0$ and a single imported case dies out.

Used in the real world· Medicine: how long to flip a proton in an MRI scanner

An MRI scanner at $1.5$ T uses protons in the water of the body. Their Larmor frequency is $63.87$ MHz, so the radio-frequency coil is tuned to it. A pulse of rotating-field amplitude $B_1 = 10\ \mu$ T tips the magnetisation by $90^\circ$ in $587\ \mu$ s (computed above), and by $180^\circ$ in $1.17$ ms. Tissue contrast comes from how quickly the tipped magnetisation relaxes back, which adds damping (the real part) to the pure rotation of the unitary evolution. A pulse that is slightly off resonance tips by less, by the Rabi formula, which is why the frequency and the field strength must be uniform across the body.

Traps

Watch out· Normal modes need a symmetric matrix

For equal masses the matrix $M$ in $x'' = -Mx$ is symmetric and the proof works. For unequal masses $m_i$ the equation is $x'' = -\mathcal M^{-1}Kx$ with $\mathcal M = \operatorname{diag}(m_i)$ and $K$ the symmetric stiffness matrix; the matrix $\mathcal M^{-1}K$ is not symmetric, its eigenvectors are not orthogonal in the dot product, and you must use the inner product weighted by the masses. The easy fix is the substitution $y = \mathcal M^{1/2}x$ , which gives $y'' = -\mathcal M^{-1/2}K\mathcal M^{-1/2}y$ with a symmetric matrix.

Watch out· The beat frequency is not a natural frequency

When one mass is pushed, mass 1 oscillates at the mean frequency $\bar\omega = 6.63$ rad/s with an envelope $\cos\Omega t$ , $\Omega = \frac{\omega_2 - \omega_1}2$ . The system has no normal mode at $\bar\omega$ . The slow envelope $\cos\Omega t$ has angular frequency $\Omega$ , but the energy goes like $\cos^2\Omega t$ , so the energy of mass 1 returns after $\frac\pi\Omega = \frac{2\pi}{\omega_2 - \omega_1}$ , twice as long as the first transfer takes. Use the eigenvalues for the natural frequencies, and sum-to-product for the beats.

Watch out· Current peaks at ω0, charge peaks a little lower

The current amplitude $\frac{E_0}{\lvert Z\rvert}$ is largest exactly at $\omega_0 = \frac{1}{\sqrt{LC}}$ . The charge on the capacitor (and so the voltage across it) has amplitude $\frac{E_0}{L\sqrt{(\omega_0^2 - \omega^2)^2 + 4\alpha^2\omega^2}}$ , which peaks at $\omega = \sqrt{\omega_0^2 - 2\alpha^2}$ . For a high- $Q$ circuit the difference is negligible ( $999.984$ kHz against $1000$ kHz in the tuner), but for $Q$ below $\frac1{\sqrt2} \approx 0.71$ (that is, $\omega_0^2 < 2\alpha^2$ ) there is no peak in the charge at all.

Watch out· R0 above 1 means growth, not the size of the epidemic, and not the herd immunity at the end

$R_0$ is a ratio. The speed of the start is $\gamma(R_0 - 1)$ , which depends on $\gamma$ as well. The fraction who are ever infected is $1 - s_\infty$ with $s_\infty$ from the final-size equation, and it is much larger than the herd immunity fraction $1 - \frac1{R_0}$ : for $R_0 = 2$ the herd immunity threshold is $50\%$ but $79.7\%$ are infected, because the epidemic coasts past the point where the infection stops growing. And $R_0$ describes a fully susceptible population; with a fraction $p$ immune the relevant number is $R_0(1 - p)$ .

Watch out· Quantum evolution: the exponent is a matrix, with the factor i, and the angle is half

$e^{-iHt/\hbar}$ is NOT the matrix of $e^{-iH_{jk}t/\hbar}$ (the trap of 7.4 again). For $H = \frac\hbar2N$ with $N^2 = \Omega_R^2I$ the angle in the cosine and sine is $\frac{\Omega_Rt}{2}$ , not $\Omega_Rt$ ; this is why the probability $\sin^2\frac{\Omega_Rt}2$ oscillates at the angular frequency $\Omega_R$ and a $\pi$ -pulse has $\Omega t = \pi$ . Check any answer for $\lvert\psi_1\rvert^2 + \lvert\psi_2\rvert^2 = 1$ .

Problems

checkcompare the worked example

Two masses of $1$ kg lie between walls with springs $k = 4$ N/m to the walls and a coupling spring $k_c = 1$ N/m. Find the normal modes and the frequencies.

Solution

$M = \begin{bmatrix}5 & -1\\-1 & 5\end{bmatrix}$ . The vector $(1, 1)$ gives $M(1, 1) = (4, 4)$ : eigenvalue $4$ , so $\omega_1 = 2$ rad/s (in step; the coupling spring is not stretched, so this is $\sqrt{k/m}$ ). The vector $(1, -1)$ gives $(6, -6)$ : eigenvalue $6$ , $\omega_2 = \sqrt6 = 2.449$ rad/s (in opposition; $\sqrt{(k + 2k_c)/m}$ ). Check: trace $10 = 4 + 6$ , determinant $24 = 4\cdot6$ .

check7.2

An RLC circuit has $L = 0.1$ H and $C = 100\ \mu$ F. Classify the free motion for $R = 100\ \Omega$ , and find the decay rates.

Solution

$2\sqrt{L/C} = 63.2\ \Omega < R = 100\ \Omega$ : overdamped. With $\alpha = \frac{R}{2L} = 500$ and $\omega_0^2 = \frac{1}{LC} = 10^5$ : $r = -500 \pm\sqrt{250\,000 - 100\,000} = -500 \pm 387.3$ , that is $-112.7$ and $-887.3$ per second. The charge is a sum of two decaying exponentials, and the slower one (time constant $\frac{1}{112.7} = 8.9$ ms) dominates after the first few milliseconds. (For $R = 20\ \Omega$ the same circuit is underdamped with $\omega_d = 300$ rad/s.)

check7.6

An SIR epidemic has $\beta = 0.5$ and $\gamma = 0.25$ per day. Find $R_0$ and the herd immunity threshold. Does an epidemic start if $60\%$ of the population is immune? At what rate does it grow when $20\%$ is immune, and what is the doubling time?

Solution

$R_0 = \frac{0.5}{0.25} = 2$ , and the threshold is $1 - \frac12 = 50\%$ . With $60\%$ immune, $s^* = 0.4$ and $R_0s^* = 0.8 < 1$ : the growth eigenvalue $\gamma(R_0s^* - 1) = 0.25\times(-0.2) = -0.05$ per day is negative, so no epidemic. With $20\%$ immune, $s^* = 0.8$ , $R_0s^* = 1.6$ , and the rate is $0.25\times0.6 = 0.15$ per day; the doubling time is $\frac{\ln2}{0.15} = 4.6$ days.

check7.7 theorem

A qubit has $H = \frac\hbar2\Omega\sigma_x$ with $\Omega/2\pi = 5$ MHz, and starts in level 1. What is the probability of level 2 after $50$ ns?

Solution

On resonance, $P = \sin^2\frac{\Omega t}{2}$ . $\Omega = 2\pi\times5\times10^6 = 3.1416\times10^7$ rad/s and $\Omega t = 3.1416\times10^7\times5\times10^{-8} = 1.5708 = \frac\pi2$ . So $P = \sin^2\frac\pi4 = \frac12$ . The pulse is a $\frac\pi2$ -pulse.

use3.1 and the chain formula

Three masses of $1$ kg lie between walls with all springs $k = 9$ N/m. Give the frequencies and the motion when $x(0) = (1, 0, -1)$ and the masses start at rest.

Solution

Frequencies $\omega_j = 2\sqrt{k/m}\sin\frac{j\pi}{8} = 6\sin\frac{j\pi}{8}$ : $2.296$ , $4.243$ , $5.543$ rad/s. The vector $(1, 0, -1)$ is the second mode shape (the middle mass stays still and the outer ones oppose). So the motion is the pure mode $x(t) = (1, 0, -1)\cos(4.243t)$ , with $\omega_2 = 6\sin\frac\pi4 = 3\sqrt2 = 4.243$ and period $1.48$ s. Check: $M(1, 0, -1) = 9(2, 0, -2) = 2\cdot9\,(1, 0, -1)$ , so $\omega^2 = 18$ and $\omega = \sqrt{18} = 4.243$ .

usethe impedance formula

A series RLC circuit has $L = 0.5$ H, $C = 2\ \mu$ F, $R = 25\ \Omega$ and is driven by $E_0 = 10$ V. Find $\omega_0$ , $Q$ , the bandwidth, the current at resonance, and the current at $\omega = 1050$ rad/s.

Solution

$\omega_0 = \frac{1}{\sqrt{0.5\times2\times10^{-6}}} = 1000$ rad/s. $Q = \frac{\omega_0L}{R} = \frac{500}{25} = 20$ . Bandwidth $\Delta\omega = \frac RL = 50$ rad/s. At resonance $\lvert Z\rvert = R$ , so the current amplitude is $\frac{10}{25} = 0.40$ A. At $1050$ rad/s the reactance is $\omega L - \frac1{\omega C} = 525 - 476.2 = 48.8\ \Omega$ , so $\lvert Z\rvert = \sqrt{25^2 + 48.8^2} = 54.8\ \Omega$ and the amplitude is $\frac{10}{54.8} = 0.182$ A: less than half the resonant value (the frequency is $50$ rad/s above resonance, one bandwidth).

usefinal-size equation

For an SIR epidemic with $R_0 = 2$ started from a nearly fully susceptible population, find the final fraction never infected $s_\infty$ by iterating $s \mapsto e^{-2(1 - s)}$ , and the fraction ever infected.

Solution

Start with $s = 0.5$ : the iterates are $0.368, 0.283, 0.238, 0.218, 0.209, 0.206, 0.204, 0.2036, \dots$ converging to $s_\infty = 0.2032$ (the iteration converges because the derivative at the root, $2s_\infty = 0.41$ , is below $1$ ). So $20.3\%$ are never infected and $79.7\%$ are infected at some time. Check in the equation: $\ln\frac{1}{0.2032} = 1.594$ and $2(1 - 0.2032) = 1.594$ .

useRabi formula

A qubit with Rabi frequency $\Omega/2\pi = 1$ MHz is driven with detuning $\delta/2\pi = 2$ MHz. What is the largest probability of flipping, and when is it first reached?

Solution

$\Omega_R = \sqrt{\delta^2 + \Omega^2} = 2\pi\sqrt{5}$ MHz $= 2\pi\times2.236$ MHz. $P_{\max} = \frac{\Omega^2}{\Omega_R^2} = \frac15 = 0.2$ . It is first reached when $\frac{\Omega_Rt}{2} = \frac\pi2$ , that is $t = \frac{\pi}{\Omega_R} = \frac{1}{2\times2.236\ \text{MHz}} = 0.2236\ \mu$ s.

prove3.1

Prove that for $x'' = -Mx$ with $M$ symmetric, the energy $E = \frac12\lVert x'\rVert^2 + \frac12x^{\mathsf T}Mx$ is constant along solutions.

Solution

Plan. Differentiate, and use symmetry of $M$ .

$\frac{d}{dt}\frac12x'\cdot x' = x'\cdot x'' = -x'\cdot Mx$ . And $\frac{d}{dt}\frac12x^{\mathsf T}Mx = \frac12(x')^{\mathsf T}Mx + \frac12x^{\mathsf T}Mx' = x'\cdot Mx$ , because $x^{\mathsf T}Mx' = (Mx)\cdot x'$ by symmetry of $M$ ( $x^{\mathsf T}Mx' = (x')^{\mathsf T}M^{\mathsf T}x = (x')^{\mathsf T}Mx$ ). The two derivatives cancel, so $\dot E = 0$ . $\blacksquare$

prove3.3 and 7.4

Prove that if $H$ is Hermitian then $e^{-iHt/\hbar}$ is unitary for every real $t$ , and deduce that $\lVert\psi(t)\rVert$ is constant.

Solution

Plan. Take the conjugate transpose of the exponential series.

Write $A = -\frac{i}{\hbar}tH$ . Then $A^* = \frac{i}{\hbar}tH^* = \frac i\hbar tH = -A$ . The adjoint of a power is the power of the adjoint, and the adjoint is continuous, so $(e^A)^* = \sum\frac{(A^*)^k}{k!} = e^{A^*} = e^{-A}$ . Since $A$ and $-A$ commute, $e^Ae^{-A} = e^0 = I$ (7.4). So $(e^A)^*e^A = e^{-A}e^A = I$ : $e^A$ is unitary. Then $\lVert e^A\psi\rVert^2 = \psi^*(e^A)^*e^A\psi = \psi^*\psi = \lVert\psi\rVert^2$ . $\blacksquare$

proveKermack and McKendrick 1927; new

Prove that if $R_0s_0 \le 1$ then the number infected $i(t)$ decreases for all $t \ge 0$ in the SIR model.

Solution

Plan. Show that $i' \le 0$ at all times by comparing $s$ with $s_0$ .

From $s' = -\beta si \le 0$ , $s$ is non-increasing, so $s(t) \le s_0$ . Then $i' = \beta si - \gamma i = \gamma i(R_0s - 1) \le \gamma i(R_0s_0 - 1) \le 0$ , because $i \ge 0$ . So $i$ is non-increasing. (If $R_0s_0 < 1$ it decreases strictly, and the infection dies out.) $\blacksquare$

prove7.2 and the impedance theorem

Prove that the half-power bandwidth of a series RLC circuit, defined by $\lvert Z\rvert^2 = 2R^2$ , is exactly $\Delta\omega = \frac RL$ .

Solution

Plan. Solve $\lvert\omega L - \frac1{\omega C}\rvert = R$ for $\omega$ and subtract the two positive roots.

$\lvert Z\rvert^2 = R^2 + X^2$ with $X = \omega L - \frac{1}{\omega C}$ , so $\lvert Z\rvert^2 = 2R^2$ means $X = \pm R$ . Multiply $\omega L - \frac{1}{\omega C} = \pm R$ by $\omega C$ : $LC\omega^2 \mp RC\omega - 1 = 0$ , so $\omega = \frac{\pm RC + \sqrt{R^2C^2 + 4LC}}{2LC} = \pm\frac{R}{2L} + \sqrt{\frac{R^2}{4L^2} + \frac{1}{LC}}$ . The positive root for the upper sign minus the positive root for the lower sign is $\frac R{2L} - \left(-\frac R{2L}\right) = \frac{R}{L}$ . $\blacksquare$

stretchMurray 10

(The SEIR model.) A latent period is added: infected people are first ‘exposed’ ( $e$ ), then infectious ( $i$ ): $s' = -\beta si$ , $e' = \beta si - \sigma e$ , $i' = \sigma e - \gamma i$ . Linearize at the disease-free state $s = 1$ and show that an epidemic starts if and only if $R_0 = \beta/\gamma > 1$ , the same as before. Compute the early growth rate for $\beta = 0.5$ , $\sigma = 0.25$ , $\gamma = 0.2$ per day, and compare with the SIR growth rate.

Solution

Plan. Linearize the $(e, i)$ part and find when it has a positive eigenvalue.

With $s \approx 1$ the equations for the exposed and infectious are $e' = \beta i - \sigma e$ , $i' = \sigma e - \gamma i$ , that is $\begin{bmatrix}e\\i\end{bmatrix}' = \begin{bmatrix}-\sigma & \beta\\\sigma & -\gamma\end{bmatrix}\begin{bmatrix}e\\i\end{bmatrix}$ . The trace is $-(\sigma + \gamma) < 0$ and the determinant is $\sigma\gamma - \sigma\beta = \sigma(\gamma - \beta)$ . By the table of 7.6 the matrix has an eigenvalue with positive real part exactly when its determinant is negative, that is $\beta > \gamma$ : $R_0 > 1$ . The largest eigenvalue is $\lambda = \frac{-(\sigma + \gamma) + \sqrt{(\sigma - \gamma)^2 + 4\sigma\beta}}{2}$ . For $\beta = 0.5$ , $\sigma = 0.25$ , $\gamma = 0.2$ : $R_0 = 2.5$ , and $\lambda = \frac{-0.45 + \sqrt{0.0025 + 0.5}}{2} = \frac{-0.45 + 0.7089}{2} = 0.129$ per day (doubling time $5.4$ days). The SIR model with the same $\beta$ and $\gamma$ grows at $\beta - \gamma = 0.3$ per day (doubling in $2.3$ days). The latent period leaves the threshold unchanged and slows the early growth. $\blacksquare$

stretch2.8 and a trigonometric identity

Prove that for $N$ masses in a row between two walls with all springs $k$ and all masses $m$ , the vector $v^{(j)} = \big(\sin\frac{ji\pi}{N+1}\big)_{i=1}^N$ is an eigenvector of $M = \frac km\,\text{tridiag}(-1, 2, -1)$ with eigenvalue $\frac{2k}{m}\big(1 - \cos\frac{j\pi}{N+1}\big) = \frac{4k}{m}\sin^2\frac{j\pi}{2(N+1)}$ .

Solution

Plan. Apply $M$ to the vector, use the sum formula for sine, and watch the boundary terms vanish.

Put $a = \frac{j\pi}{N+1}$ and $v_i = \sin(ia)$ , $i = 0, 1, \dots, N+1$ , so that $v_0 = \sin0 = 0$ and $v_{N+1} = \sin(j\pi) = 0$ : these two ‘virtual’ entries are exactly the wall conditions. For $1 \le i \le N$ : $(Mv)_i = \frac km(2v_i - v_{i-1} - v_{i+1})$ (with $v_0 = v_{N+1} = 0$ the first and last rows are correct too). By $\sin(ia - a) + \sin(ia + a) = 2\sin(ia)\cos a$ , we get $v_{i-1} + v_{i+1} = 2\cos a\,v_i$ . So $(Mv)_i = \frac km(2 - 2\cos a)v_i$ , which is the eigenvalue equation with $\lambda = \frac{2k}{m}(1 - \cos a)$ . Since $1 - \cos a = 2\sin^2\frac a2$ , $\lambda = \frac{4k}{m}\sin^2\frac{j\pi}{2(N+1)}$ and $\omega_j = 2\sqrt{k/m}\sin\frac{j\pi}{2(N+1)}$ . The vectors are non-zero for $j = 1, \dots, N$ and have distinct eigenvalues, so they are all the modes. $\blacksquare$

Summary

Coupled springs: $x'' = -Mx$ with $M$ symmetric positive definite. The eigenvectors are the normal modes (patterns that oscillate at one frequency), $\omega_j = \sqrt{\text{eigenvalue}}$ , and every motion is a sum of modes (beats appear when two frequencies are close). The energy of each mode is constant.
An RLC circuit obeys $Lq'' + Rq' + q/C = E$ : the oscillator equation with $L \leftrightarrow m$ , $R \leftrightarrow c$ , $1/C \leftrightarrow k$ . It is always asymptotically stable ( $\tau = -R/L$ , $\Delta = 1/LC$ ), underdamped if $R < 2\sqrt{L/C}$ , and $\dot W = IE - RI^2$ for its energy.
Driven by $E_0\cos\omega t$ the steady current is $\frac{E_0}{\lvert Z\rvert}\cos(\omega t - \varphi)$ with the impedance $Z = R + i(\omega L - \frac{1}{\omega C})$ ; resonance at $\omega_0 = 1/\sqrt{LC}$ , bandwidth $R/L$ , quality factor $Q = \omega_0L/R$ .
The SIR model $s' = -\beta si$ , $i' = \beta si - \gamma i$ has a line of disease-free equilibria. The Jacobian there has eigenvalues $0$ and $\gamma(R_0s^* - 1)$ , so an epidemic starts exactly when $R_0s^* > 1$ ; with $R_0 = \beta/\gamma$ and herd immunity at $1 - 1/R_0$ .
The quantity $i + s - \frac{1}{R_0}\ln s$ is conserved, giving the peak $1 - \frac{1 + \ln R_0}{R_0}$ and the final size (solve $\ln\frac{s_0}{s_\infty} = R_0(1 - s_\infty)$ ): for $R_0 = 3$ , $30\%$ at the peak and $94\%$ ever infected.
A two-level quantum system evolves by $\psi(t) = e^{-iHt/\hbar}\psi(0)$ . For Hermitian $H$ the matrix is unitary, so probability is conserved, and for $H = \frac\hbar2(\delta\sigma_z + \Omega\sigma_x)$ the Rabi formula gives $P(t) = \frac{\Omega^2}{\Omega_R^2}\sin^2\frac{\Omega_Rt}{2}$ : a $\pi$ -pulse at $\Omega/2\pi = 10$ MHz lasts $50$ ns.
In practice: turn the law into $x' = Ax$ (or $x'' = -Mx$ ), find the eigenvalues, and read their meaning: imaginary part, frequency; real part, rate; sign of one eigenvalue, a threshold; and keep an energy or a conserved quantity as a check.

Watch and read

The same ideas explained by someone else. Use these when an explanation here does not land.

7.7 Applications to springs, circuits, epidemics and qubits

The picture

The idea, in words

Definitions and theorems

Coupled springs and normal modes

Normal mode

Normal-mode decomposition

RLC circuits

RLC circuit

A passive circuit loses energy

Steady-state response of the driven RLC circuit

Quality factor

The SIR epidemic

SIR model

Basic reproduction number

Disease-free equilibrium

The epidemic threshold

Peak and final size of an epidemic

A two-level quantum system

Hamiltonian

Two-level system

Unitary evolution

The Rabi formula

Where it is used

Traps

Problems

Summary

Watch and read

Where this sits

Comes from

Leads to

Same idea elsewhere