Stabilize high-order Maxwellian initialization

[rogeriojorge]

Summary

• generate drifting-Maxwellian Hermite coefficients with their adjacent-mode recurrence
• avoid separately evaluating exponentially large powers and factorials
• retain the existing Fourier transform, public API, shapes, and low-order values

Root cause

compute_C_nmp formed 2 ** n, factorial(n), alpha ** (n + 1), and
(U - u) ** n independently. At high Hermite order these intermediates
overflow even when their combined coefficient is small, so inf / inf
produces NaNs in the initial state before integration. An Orszag--Tang state
with Nn=112 reproduced the failure.

The replacement evaluates the mathematically equivalent one-dimensional
recurrence

c_0 = 1 / alpha
c_(n+1) = c_n * ((U - u) / alpha) * sqrt(2 / (n + 1))

for each velocity axis and multiplies the three broadcast coefficient tables.
This avoids large cancelling intermediates and reduces the amount of repeated
power/factorial work.

Validation

• pytest -q: 4 passed
• fatal flake8 checks: 0
• recurrence agrees with the previous closed form at low order to 1e-13
• Nn=256, Nm=4, Np=4 initialization is finite
• 129x129, 112x4x4 Orszag--Tang initialization changed from non-finite to finite
• on that case, cold/warm initialization was 0.672/0.587 s versus
  0.679/0.606 s for the previous non-finite expression

This PR is independent of the parallelization stack and targets main.

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.

Files with missing lines	Coverage Δ
spectrax/_initialize_maxwellian.py	100.00% <100.00%> (+62.96%)	⬆️

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.