Unwanted sine modes appearing in convolution using FFT

[pw0908]

Hi All, I don't know if this is an issue with FastTransforms.jl specifically or something that I've done, but any guidance that could be provide will be appreciated.

To give a little bit of background, I'm trying to perform the following convolution (nV, r and r' are vectors):

$$n_\mathbf{V}(r) = \int d\mathbf{r}' \rho (r') \delta(|r-r'|-R) \frac{r-r'}{|r-r'|}$$

Where rho(r) is a profile which looks like (I've made sure it's periodic):

I have obtained an analytical expression for the fourier transform of the last two terms in the integral as:

$$\hat{w} = 4\pi i \frac{k}{|k|^3} (\sin(|k|R)-|k|R\cos(|k|R))$$

If I use just regular Float64, the resulting nV profile looks fine:

Except when I zoom in:

While these values are small, in a later part of my code, I need to divide this profile by a small number which makes this 'noise' blow up and affects calculations downstream.

I tried using FastTransforms.jl since it would let me go to higher precision. However, even when I use BigFloat, although it removes the noise, some unphysical oscillations remain:

Im no longer certain that these remaining oscillations are due to floating point errors. I also call these oscillations unphysical as, aside from the nature of this convolution integral, if I perform this convolution integral manually, the oscillations disappear completely (and I'm able to obtain the integral to higher precision).

Am I doing something wrong in the way I'm performing the convolution integral? Thanks!

I've attached an MWE below:

using FFTW, FastTransforms

# Generating the density profile
tanh_prof(x,start,stop,shift,coef) = 1/2*(start-stop)*tanh((x-shift)*coef)+1/2*(start+stop)

ub = 10
lb = -10
z = LinRange(lb,ub,4000)
rho1 = 1e3
rho2 = 1e-12
rho = @. tanh_prof(z,rho1,rho2,(ub/4+3*lb/4),20)*(z<=0) +
         tanh_prof(z,rho2,rho1,(3*ub/4+lb/4),20)*(z>0)

# Obtaining the fourier transform of w_hat
freq = fftfreq(4000,4000/20)
R = 1/2*2pi

weight_hat = @. 0.0 - 4π*im*freq/abs(freq)^3*(sin(abs(freq)*R)-R*abs(freq)*cos(abs(freq)*R)) *(freq != 0.0)

# Performing the convolution integral
rho_hat = fft(rho)
I_hat = @. rho_hat*weight_hat
nV = real(ifft(I_hat))

[pw0908]

One thing to add: increasing the number of points does not remove these oscillations and doesn't change their amplitude.