Potential issue with `_make_hermitian()` causing non-real IFFT output

[karthikviswanathn]

Hi @steven-murray, thank you for maintaining this library. I find it really helpful for my research! I came across a potential issue in the code when I noticed that the Fourier transform of pb.delta_k() had giving imaginary parts. Since pb.delta_k()is intended to be Hermitian, the fourier transform must be real, which did not seem to be the case. Can you let me know if this is intended behaviour or a bug that should be fixed?

Summary

gauss_hermitian() may not produce perfectly Hermitian arrays, leading to IFFT outputs with imaginary components.

Expected vs Observed

• Expected: For Hermitian FFT arrays, IFFT should produce real output (within numerical precision)
• Observed: IFFT(gauss_hermitian()) has max imaginary component ~0.015

Reproducible Test

import numpy as np
import powerbox as pbox

pb = pbox.PowerBox(N=32, dim=2, pk=lambda k: k**(-2.0),
                    boxlength=32.0, seed=42, ensure_physical=False)
gh = pb.gauss_hermitian()

# Check Hermitian property: gh[i,j] = conj(gh[(-i)%N, (-j)%N])
N = gh.shape[0]
max_err = max(abs(gh[i,j] - np.conj(gh[(-i)%N, (-j)%N]))
            for i in range(N) for j in range(N))

print(f"Hermitian error: {max_err:.2e}")  # Output: 3.48e+00
print(f"Max |imag(IFFT)|: {abs(np.imag(np.fft.ifft2(gh))).max():.2e}")  # Output: 1.54e-02

Potential Cause

_make_hermitian() uses [::-1] which maps [i,j] → [N-1-i, N-1-j], but FFT Hermitian symmetry needs [i,j] → [(-i)%N, (-j)%N]. These differ at boundaries (e.g., i=0: N-1 vs 0).

Environment: PowerBox 0.8.2, NumPy 2.3.1, Python 3.11.7