Wave packets

A perfectly sinusoidal wave extends forever in space and time. It has a precise wave number and frequency, but it cannot represent a localized pulse or finite message. A wave packet combines a range of wave numbers with carefully related phases so that they add constructively in one region and mostly cancel elsewhere.

Localization requires multiple components

A one-dimensional linear packet can be written as a continuous superposition

\[\psi(x,t)=\frac{1}{2\pi}\int_{-\infty}^{\infty}\widetilde{\psi}(k,0)e^{i[kx-\omega(k)t]}\,dk\]

The complex coefficient \(\widetilde{\psi}(k,0)\) sets the amplitude and phase of wave-number component \(k\). The dispersion relation

\(\omega(k)\)

sets how that component's phase develops.

At a location where many phases line up, the sum is large. Far away, components point in different complex directions and largely cancel. The packet is therefore not a separate particle-like ingredient added to sinusoidal waves; it is the interference pattern made by them.

Two nearby waves reveal an envelope

Add two equal waves with nearby wave numbers and angular frequencies:

\[\cos(k_1x-\omega_1t)+\cos(k_2x-\omega_2t)=2\cos\left(\frac{\Delta k\,x-\Delta\omega\,t}{2}\right)\cos\left(\bar{k}x-\bar{\omega}t\right)\]

Here \(\Delta k=k_2-k_1\), \(\Delta\omega=\omega_2-\omega_1\), \(\bar{k}=(k_1+k_2)/2\), and \(\bar{\omega}=(\omega_1+\omega_2)/2\). The rapidly varying second cosine is the carrier. The slowly varying first cosine is its envelope.

Holding the envelope phase constant gives its speed:

\[v_{\mathrm{envelope}}=\frac{\Delta\omega}{\Delta k}\]

As the two components become arbitrarily close in wave number, this becomes group velocity:

\[v_g=\left.\frac{d\omega}{dk}\right|_{k_0}\]

The carrier crests move with phase velocity

\[v_p=\frac{\omega(k_0)}{k_0}\]

These velocities are equal only for special dispersion relations. A crest can move through an envelope without carrying the packet center along with it.

A Gaussian packet with defined widths

Width statements are meaningless unless they say whether width belongs to amplitude, intensity, energy density, or probability density. One convenient normalized complex packet is

\[\psi(x,0)=\frac{1}{(2\pi\sigma_x^2)^{1/4}}\exp\left[-\frac{(x-x_0)^2}{4\sigma_x^2}\right]e^{ik_0(x-x_0)}\]

With this definition, \(|\psi|^2\) is a Gaussian whose standard deviation in position is \(\sigma_x\). Using the spatial transform convention in the previous lesson, its spectral squared magnitude is Gaussian with standard deviation \(\sigma_k\) and

\[\sigma_x\sigma_k=\frac12\]

The central wave number \(k_0\) controls the carrier spacing. The spectral width controls localization. Making \(\sigma_k\) four times larger makes \(\sigma_x\) four times smaller under this definition.

The real part may be drawn as a Gaussian envelope multiplying a cosine, but calculations should retain the complex packet or both quadratures. Taking only the real part changes some normalization and spectral bookkeeping because it introduces positive and negative wave-number components.

Spectral time evolution

In a linear, time-independent medium, evolve every spectral component by a unit-magnitude phase:

\[\widetilde{\psi}(k,t)=\widetilde{\psi}(k,0)e^{-i\omega(k)t}\]

For a lossless model this multiplication preserves \(|\widetilde{\psi}(k,t)|\). The inverse transform then builds the new position-space packet.

import numpy as np

N, length = 1024, 100.0
x = (np.arange(N) - N / 2) * length / N
dx = length / N
wave_numbers = 2 * np.pi * np.fft.fftfreq(N, d=dx)
initial = np.exp(-(x**2) / 8) * np.exp(2j * x)
spectrum = np.fft.fft(initial)
frequency = wave_numbers**2 / 2  # Scaled dispersive model.
packet = np.fft.ifft(spectrum * np.exp(-1j * frequency * 2.0))
print("Peak location:", x[np.argmax(abs(packet))])

This method keeps material physics in \(\omega(k)\) and numerical representation in the grid and transform. A loss model would also change spectral magnitudes and must be stated separately.

Translation in a nondispersive medium

For a right-moving nondispersive relation \(\omega(k)=vk\), the phase factor combines with \(e^{ikx}\) to give \(e^{ik(x-vt)}\). Therefore

\[\psi(x,t)=\psi(x-vt,0)\]

The entire ideal packet translates at speed \(v\) without changing shape. A numerical calculation may still broaden it because of grid dispersion or a dissipative interpolation method; that is numerical behavior, not physical dispersion.

Dispersion and packet spreading

Near the central wave number, expand the dispersion relation:

\[\omega(k)\approx\omega_0+\omega'_0(k-k_0)+\frac12\omega''_0(k-k_0)^2+\cdots\]

The constant term advances the overall phase. The first derivative \(\omega'_0\) moves the envelope at group velocity. The second derivative \(\omega''_0\) makes different parts of the spectrum accumulate nonlinear relative phase, usually broadening and chirping the packet.

Group velocity is a local, narrow-band approximation. In strongly dispersive, absorptive, or reshaping media, the position of a packet maximum can behave strangely and need not equal the speed of information or energy. A true signal front is constrained by the full causal medium response, not by one derivative evaluated at \(k_0\).

Measuring center and width

For a nonnegative weight \(w(x,t)=|\psi(x,t)|^2\), define

\[\langle x\rangle=\frac{\int xw\,dx}{\int w\,dx},\qquad \sigma_x^2=\frac{\int(x-\langle x\rangle)^2w\,dx}{\int w\,dx}\]

For a quantum wavefunction, normalized \(|\psi|^2\) has a probability interpretation. For a classical displacement field, squared displacement is not generally the complete energy density. Use the actual kinetic and elastic energy when the desired center is an energy centroid.

On a periodic numerical domain, ordinary coordinate averages fail when a packet straddles the wrap point. Use circular statistics or unwrap the packet around a tracked center.

Packet, pulse, and envelope are related but distinct

A pulse is a disturbance localized in time, space, or both. A packet is a pulse described through a band of Fourier components. The envelope is a smooth curve used to characterize the packet's slowly varying amplitude. It is not normally another physical field added to the carrier. Multiplying a carrier by an envelope in position space corresponds to convolving their spectra, which is why stronger localization requires broader spectral content.

An envelope is clearest when its bandwidth is small compared with its central wave number. Then several carrier oscillations fit beneath it and the split into “fast carrier” and “slow envelope” is useful. For a pulse containing only one or two oscillations, different envelope definitions can give noticeably different centers and widths. Report the definition used instead of treating the drawn outline as unique.

Classical and quantum packets share Fourier mathematics but not every interpretation. A classical string displacement may be positive or negative, and its mechanical energy requires displacement derivatives and velocity. A quantum wavefunction is complex, and its squared magnitude is a probability density after normalization. A classical packet does not become a particle because it is localized, and a quantum packet is not a tiny classical wave blob following one definite path.

Measurements also have finite bandwidth. A sensor records its own filtered version of the field, so a measured packet width can combine source width, medium dispersion, detector response, and sampling. Calibrate or model those stages before attributing all broadening to propagation.

Grid, domain, and boundary constraints

With spatial spacing \(\Delta x\), the largest distinct wave-number magnitude is associated with the spatial Nyquist limit, approximately

\[|k|\le\frac{\pi}{\Delta x}\]

A carrier too close to this limit has too few samples per wavelength and is highly sensitive to numerical dispersion. The domain must also be several packet widths across. Because a DFT implies periodic boundaries, a packet leaving one side reappears on the other. That wraparound is correct for a ring and wrong for an intended open line.

initial_norm = np.sum(abs(initial) ** 2) * dx
assert 2.0 + 3 * 0.5 < 0.8 * np.pi / dx
for time in [0.0, 1.0, 2.0]:
    packet = np.fft.ifft(spectrum * np.exp(-1j * frequency * time))
    norm = np.sum(abs(packet) ** 2) * dx
    centre = np.sum(x * abs(packet) ** 2) * dx / norm
    assert np.isclose(norm, initial_norm)
    print(time, centre, norm)

Verification

Construct two Gaussian spectra with the same \(k_0\) but spectral widths in a ratio of four. Their initial position widths should have the inverse ratio. Multiply the spectrum by \(e^{-ikx_s}\) and verify that the reconstructed packet shifts by \(x_s\) without changing width or norm.

Under \(\omega=vk\), the measured center should move linearly with slope \(v\) and width should remain fixed up to numerical error. Under a relation with nonzero curvature, verify the predicted initial group velocity and observe convergence of width evolution as spatial resolution and domain size increase.

Real wave packets can experience loss, nonlinear self-interaction, scattering, boundaries, and multiple dimensions. A quantum packet additionally requires a Hamiltonian, normalization, and measurement interpretation. The Fourier construction supplies the kinematics of localization; it does not by itself choose the physical medium.

A worked Python notebook

Follow the standard deviation of a freely spreading minimum-uncertainty Gaussian with initial width one. The question for this calculation is precise: how does Gaussian position width change when we change time, while keeping the other choices fixed? The earlier sections explain the physical relationship. Here we turn one part of that relationship into a calculation that you can run, inspect, and change in Google Colab.

The example uses mass and reduced Planck constant equal to one and an effectively unbounded domain. Keep that sentence beside the results. A formula can return a number outside its useful range, and a plotting library will usually draw that number without asking whether the physical assumptions still hold. The selected interval is from 0 to 5 scaled time units. It gives us a specific domain to discuss instead of leaving the graph's range to a software default.

Use the notebook link above the lesson, run the setup cell, and then run the cells in order. The code below is also included there. NumPy handles the arrays of numerical values, SymPy keeps the equation symbolic, and Matplotlib draws the results. If these libraries or notebook cells are unfamiliar, start with Your first physics notebook.

State the relationship before choosing a graph

In this small example, q stands for time, expressed in scaled time units. The output is Gaussian position width, expressed in scaled length units. The equation below uses the numerical values in those units. Any fixed coefficients therefore carry the units needed to make the relation consistent. A value supplied in another unit must be converted before it enters this expression.

\[f(q) = \sqrt{\frac{q^{2}}{4} + 1}\]

This is a function of one changing input. It is not a statement that the physical system has only one variable. We deliberately hold the other inputs fixed so that one cause of variation is visible. Changing several inputs at once would make it harder to tell which change produced the output difference. Write down the held-fixed choices before trying a different experiment.

At 2 scaled time units, the expression gives 1.41421 scaled length units. This is our reference calculation. It is a prediction of the stated example, not a new measurement. Its displayed digits help us compare implementations. They should not be read as a claim that physical inputs are known to six significant figures.

import numpy as np
import sympy as sp
import matplotlib.pyplot as plt

q = sp.Symbol("q", real=True)
expression = sp.sympify("sqrt(1+q**2/4)", locals={"q": q})
reference_input = 2.0
reference_output = float(expression.subs(q, reference_input))
print("Input:", reference_input, "scaled time units")
print("Predicted Gaussian position width:", reference_output, "scaled length units")

sp.Symbol declares a mathematical symbol. It does not assign a numerical measurement to that symbol. subs then replaces the symbol with our chosen input, and float converts the resulting symbolic number into a standard Python numerical value. Keeping these steps separate lets us inspect the equation before numerical substitution. It also prevents a later plot from becoming the only surviving record of the model.

Compare a few cases you can explain

Before drawing hundreds of points, calculate three. The lower endpoint gives 1 scaled length units. The midpoint, at 2.5 scaled time units, gives 1.60078 scaled length units. The upper endpoint gives 2.69258 scaled length units. Read these together with the assumptions: they describe different choices of time within the same example, with the other inputs held fixed.

These three values are a useful first check, but they do not establish the shape between them. Two endpoints can be equal even when the curve bends or oscillates in between. Likewise, a large change between endpoints does not tell us whether the relationship is linear. We need the equation and intermediate samples to answer those questions.

evaluate = sp.lambdify(q, expression, modules="numpy")
inputs = np.array([0.0, 2.5, 5.0])
outputs = np.asarray(evaluate(inputs), dtype=float)
outputs = np.broadcast_to(outputs, inputs.shape)
for input_value, output_value in zip(inputs, outputs):
    print(
        f"{input_value:.6g} scaled time units -> "
        f"{output_value:.6g} scaled length units"
    )
assert np.all(np.isfinite(outputs))

lambdify turns the symbolic expression into a numerical function. We explicitly select NumPy as its numerical backend so it can accept an array. The same expression now has two useful forms: one for algebra and one for calculating many cases. The assertion checks that the chosen inputs produce finite outputs. It would catch an infinity or an undefined numerical result, but a finite value can still be physically inappropriate.

Make a graph that answers the question

The horizontal axis shows time; the vertical axis shows Gaussian position width. Every plotted point is a prediction of the same expression. We mark the reference case separately so that a reader can connect the earlier printed value to its place on the curve. The axis labels include units because a shape without a scale cannot tell us the size of a physical effect.

input_grid = np.linspace(0.0, 5.0, 301)
output_grid = np.asarray(evaluate(input_grid), dtype=float)
output_grid = np.broadcast_to(output_grid, input_grid.shape)
fig, ax = plt.subplots(figsize=(7, 4))
ax.plot(input_grid, output_grid, label="Stated model")
ax.scatter([reference_input], [reference_output], label="Reference case", zorder=3)
ax.set_xlabel("Time (scaled time units)")
ax.set_ylabel("Gaussian position width (scaled length units)")
ax.legend()
ax.grid(alpha=0.25)
fig.tight_layout()
plt.show()

Look for a sign change, a turning point, or a region where the curve becomes steep. Connect each feature to the expression rather than to the colour of the line. If the curve has a minimum or maximum, changing the input in opposite directions may have similar effects on the output near that point. If it is nearly straight over a short interval, a local linear approximation may be useful even when the complete expression is nonlinear.

Increasing the number of samples makes the line smoother and may reveal features that a coarse grid missed. It does not improve the physical inputs or enlarge the model's validity range. For quantities that can take only discrete values, use the line as a guide to the relationship and interpret only the allowed input values. A plotting function does not enforce such restrictions for you.

Measure the local sensitivity

A derivative asks a more specific question than the whole graph: how quickly does Gaussian position width change near the reference input? Its units are the output unit divided by the input unit. For this calculation, the derivative at 2 scaled time units is 0.353553 in those ratio units. Locally, the output increases as the input increases. This statement is local; the sign or size of the slope can be different elsewhere on the curve.

\[f'(q) = \frac{q}{4 \sqrt{\frac{q^{2}}{4} + 1}}\]

For a small input change, multiply the derivative by that change to estimate the output change. Then compare that estimate with a fresh evaluation of the full expression. The distinction between these two calculations matters: the derivative gives a tangent approximation, while substitution evaluates the chosen expression at the new point.

slope_expression = sp.diff(expression, q)
slope_at_reference = float(slope_expression.subs(q, reference_input))
changed_input = 2.025
linear_prediction = reference_output + slope_at_reference * (
    changed_input - reference_input
)
full_prediction = float(evaluate(changed_input))
print("Symbolic derivative:", slope_expression)
print("Tangent prediction:", linear_prediction, "scaled length units")
print("Full expression:", full_prediction, "scaled length units")
print(
    "Approximation difference:",
    full_prediction - linear_prediction,
    "scaled length units",
)

Here the changed input is 2.025 scaled time units. The tangent estimate is 1.42305 scaled length units, while the full expression gives 1.42308 scaled length units. If the two values are close, it means the chosen input change is small enough for this local approximation at this point. It does not mean the entire model is linear. If the slope is zero, a quadratic or higher-order term may supply the leading change instead.

This is also the beginning of uncertainty propagation. If the input has a small stated standard uncertainty, the magnitude of the slope times that uncertainty estimates one contribution to output uncertainty. Do not automatically use the same rule for a wide interval, a discontinuity, a turning point with zero first derivative, or several correlated uncertain inputs. Those situations need more of the model than one slope.

Check the derivative with numerical differences

We can estimate the same slope from values on either side of the reference point. Subtract the lower output from the upper output and divide by twice the input step. Parentheses are essential: the entire output difference must be divided by the entire input interval. Without them, Python can execute a different mathematical expression while still producing a plausible number.

The calculation below uses three decreasing step sizes. Each estimate is compared with the symbolic derivative, so the reference is independent of the numerical difference formula. Near an endpoint, it is acceptable here to evaluate a small extension only because the displayed mathematical expression is smooth there; such a check does not assert that extrapolated values describe additional physical cases.

input_span = 5.0
for fraction in [1e-2, 1e-3, 1e-4]:
    step = input_span * fraction
    upper = float(evaluate(reference_input + step))
    lower = float(evaluate(reference_input - step))
    numerical_slope = (upper - lower) / (2 * step)
    difference = numerical_slope - slope_at_reference
    print(f"step={step:.6g}: slope={numerical_slope:.6g}, difference={difference:.6g}")

For a smooth function, reducing the step often reduces the approximation error at first. Very small steps can eventually make subtraction of nearly equal floating-point values less reliable. Therefore “smaller” is not an unlimited guarantee of “better.” Compare the observed behaviour over a useful range of steps instead of choosing the smallest number the computer can represent.

An absolute slope difference is useful when the expected slope is near zero. A relative difference would divide by a very small reference and could look enormous even when the absolute discrepancy is tiny. Away from zero, a relative difference can help compare quantities with very different scales. The tolerance should reflect the question and the numerical method, rather than a universal decimal threshold copied into every test.

Change the example and preserve the reasoning

Make one copy of the reference calculation before editing it. Choose a second value of time inside the stated interval, predict whether Gaussian position width will rise or fall, and then evaluate the expression. Place the prediction and the result in a text cell. If they disagree, identify whether the mistake was in the sign, the unit conversion, the algebra, or the assumption about how the curve behaves.

Next choose a shorter plotting interval around the reference point. Compare the local curve with its tangent approximation. Then widen the interval again. The difference shows why a useful local estimate need not describe the whole range. Keep the axis labels and numerical scale visible when comparing the plots; automatic axis limits can make two different effects appear equally large.

Finally, return to the physical limit stated at the start: The example uses mass and reduced Planck constant equal to one and an effectively unbounded domain. Describe one change to the setup that would violate that limit. State what new input, equation, or observation would be needed before trusting an answer in the changed situation. This last step connects the Python experiment to the physics instead of treating successful execution as the only criterion for a useful result.

Before saving, restart the runtime and run all cells in order. A result that depends on an earlier, now-deleted input cell is not reproducible from the saved notebook. Keep unrounded values for later calculations, round only the reported output, and label any new data as measured, assumed, or generated. The graph and printed values should remain traceable to the same expression and the same choices of inputs.

Recover the output change by integrating the slope

There is another way to check that the expression and its derivative describe the same relationship. Add up the small output changes predicted by the slope across the chosen input interval. In the continuous limit, this integral equals the output at the upper endpoint minus the output at the lower endpoint. The result has units of scaled length units: the input unit in the integration cancels the input unit in the denominator of the slope.

This statement is about the mathematics of this particular relationship. It does not assign a new physical meaning to an arbitrary area under a graph. For example, an area under force versus distance represents work, while an area under a different pair of quantities may represent something else or have no standard name. Always carry the axis units through the calculation before naming an integral.

slope_function = sp.lambdify(q, slope_expression, modules="numpy")
endpoint_change = float(evaluate(5.0) - evaluate(0.0))
for points in [101, 201, 401]:
    grid = np.linspace(0.0, 5.0, points)
    slopes = np.asarray(slope_function(grid), dtype=float)
    slopes = np.broadcast_to(slopes, grid.shape)
    integrated_change = np.trapezoid(slopes, grid)
    print(
        points,
        "integrated change:",
        integrated_change,
        "endpoint change:",
        endpoint_change,
    )

The trapezoidal rule joins neighbouring slope samples with straight segments and adds their signed areas. A negative contribution remains negative; replacing every slope by its magnitude would calculate a different quantity. Increasing the sample count provides a practical refinement check. If the values do not settle toward the endpoint difference, first inspect the interval, singularities, and array ordering.

Do not conclude that two agreeing calculations prove the model is physically complete. Both use the same selected relationship, so an omitted effect would be omitted from both. Their agreement is evidence that the numerical tools implement that relationship consistently. Comparing with an independent experiment or a more appropriate physical model is a separate test.

Keep this comparison in the notebook when you change the expression. It is particularly useful after adding a coefficient, changing a sign convention, or converting an input unit. A code cell that produces a smooth graph can conceal these mistakes; a second calculation of the same output change gives you another place to notice them.

Sources

Share to: