Plane Wave Scattering by a Rigid Cylinder

1. Governing equation

For time-harmonic acoustic waves at angular frequency ω, the pressure satisfies the Helmholtz equation:

∇²p + k²p = 0, k = ω/c

where c is the sound speed and k is the wavenumber. In polar coordinates (r, θ), solutions separable in r and θ take the form Z_n(kr) cos(nθ) or Z_n(kr) sin(nθ), where Z_n is any Bessel function of order n.

2. Incident wave expansion (Jacobi–Anger identity)

A plane wave propagating in the +x direction has the spatial part p_inc = e^ikx = e^{ikr cosθ}. The Jacobi–Anger identity expands this into cylindrical harmonics [1, Ch. 7]:

e^{ikr cosθ} = ∑_n=0^∞ ε_n iⁿ J_n(kr) cos(nθ)

where J_n is the Bessel function of the first kind of order n, and ε_n is the Neumann factor: ε₀ = 1, ε_n = 2 for n ≥ 1. Only cosine terms appear because the incident wave and cylinder geometry are symmetric about the x-axis.

3. Scattered field and boundary condition

The scattered field must consist of outgoing waves at infinity, represented by Hankel functions of the first kind H_n⁽¹⁾. (Note: this derivation uses the e^−iωt time convention, under which H_n⁽¹⁾, with asymptotic behavior e^+ikr/√r, represents outgoing waves satisfying the Sommerfeld radiation condition. Texts using e^+iωt would use H_n⁽²⁾ instead.)

p_scat = ∑_n=0^∞ ε_n iⁿ A_n H_n⁽¹⁾(kr) cos(nθ)

The rigid (Neumann) boundary condition ∂p/∂r = 0 at r = a requires, for each order n independently:

ε_n iⁿ [J_n′(ka) + A_n H_n^(1)′(ka)] = 0

Solving for the scattering coefficient:

A_n = −J_n′(ka) / H_n^(1)′(ka)

4. Total field

Combining the incident and scattered contributions:

p(r,θ) = ∑_n=0^N ε_n iⁿ [J_n(kr) + A_n H_n⁽¹⁾(kr)] cos(nθ)

The animated display shows Re[p(r,θ) e^−iωt] at each time step, where ωt advances continuously. The modulus display shows |p(r,θ)|, which is constant in time.

Implementation note: In the demo, the incident plane wave is evaluated directly as e^ikx (a single cosine and sine per pixel) rather than via the Jacobi–Anger series expansion. This avoids truncation artifacts in the plane wave at high ka, where many terms would be needed for the series to converge far from the cylinder. The N-term series is used only for the scattered wave, which converges rapidly at N ≈ ⌈ka⌉ + 10. The total field computed is:

p = e^ikx + ∑_n=0^N−1 ε_n iⁿ A_n H_n⁽¹⁾(kr) cos(nθ)

which is mathematically equivalent to the combined form above.

5. Bessel function derivatives

Derivatives of any cylindrical Bessel function Z_n (whether J_n, Y_n, or H_n⁽¹⁾) are computed via the recurrence [4, §9.1.27]:

Z_n′(z) = Z_n−1(z) − (n/z) Z_n(z)

For n = 0: Z₀′(z) = −Z₁(z). Since H_n⁽¹⁾ = J_n + iY_n, the derivative H_n^(1)′ = J_n′ + iY_n′, making A_n a complex number computed via complex division.

6. Numerical computation of Bessel functions

The demo computes J_n(x) and Y_n(x) for integer order n and real argument x without external libraries. Seed values J₀, J₁, Y₀, Y₁ are computed from rational polynomial approximations tabulated in Abramowitz & Stegun [4, §9.4]. Higher orders use:

• J_n(x): Miller’s backward recurrence from a high starting index, normalized using the identity J₀(x) + 2∑_k=1^∞ J_2k(x) = 1. Backward recurrence is stable because J_n is the minimal (subdominant) solution of the three-term recurrence.
• Y_n(x): Forward recurrence from Y₀, Y₁ via Y_n+1(x) = (2n/x)Y_n(x) − Y_n−1(x). Forward recurrence is stable for Y_n because it is the dominant solution.

7. Series convergence

The series converges rapidly because the scattering coefficients A_n → 0 for n ≫ ka (the numerator J_n′(ka) vanishes while the denominator grows). A safe truncation is N ≈ ⌈ka⌉ + 10 [1, Ch. 8]. The demo lets you reduce N below this to observe convergence artifacts. For real-time performance, Bessel values J_n(kr), Y_n(kr) are precomputed on a radial lookup table of 512 bins, and cos(nθ) is evaluated via Chebyshev recurrence rather than repeated trigonometric calls.