Where to apply zero-padding for convolution dealiasing and appropriate scale

After hours of browsing the DSP posts and resources online, I still struggle to understand why my code diverges when I activate zero-padding dealiasing. When I deactivate it, everything 'works well', and I pinned down the issue to the points I raise below.

I'm doing a loop where I

1. retrieve from the solution of a linear system some Fourier coefficients in frequency space - say $$\hat{u}_k$$ and $$\hat{v}_k$$ - for each $$k \in [-K,K]$$
2. apply backward DFT to go to physical space.
3. form product $$uv$$
4. apply forward DFT to go back to frequency space and broadcast the de-aliased product to update the r.h.s. of the linear system for the next iteration

The Orszag "3/2" zero-padding rule ensures dealiasing with such quadratic linearities $$uv$$. I extended accordingly my DFT input/output arrays at steps 2. and 4 (from $$N$$ to $$N_\text{pad}$$) and made sure the zeros are inserted at the right location.

• Is full dealiasing achieved at step 3? In other words, should I only apply zero-padding during my backward DFT at step 2? Should my forward DFT call be of size $$N$$ only?

• It seems that there (p18) and there (p4) they apply a factor 3/2 = $$N_\text{pad}/N$$. Can someone explain mathematically or in physical terms the reason for such a scaling? I grasp that the 'numerical' FFT is of size $$N_\text{pad}$$ whereas you're interested in a convolution of a size $$N$$ only but I'm still puzzled. In this open-source library there's no such scaling*.

*I'm aware FFTW returns an unnormalized DFT in contrast to MATLAB or numpy where ifft is normalized.