# How to set up the right model for solving diffraction numerically?

The diffraction of waves in free space is governed by the sourceless Helmholtz equation $$\nabla^2u +k^2u=0$$ There is a well stablished solution called Kirchhoff´s diffraction integral (KDI) and it´s deduction is straightforward from greens theorem, the solution is $$u(r) = \int_{S} \left [ u \space \partial_n \left(\frac{e^{iks}}{s} \right) - \left(\frac{e^{iks}}{s} \right) \space \partial_n u\space \right] dS$$

Where $$S$$ is a close surface containing the aperture and an imaginary boundary at infinity that fulfills Sommerfeld radiation condition and $$s$$ is the distance between any point in the surface and the point $$r$$ inside the enclosed volume.

Now I am trying to solve the Helmholtz equation numerically using the Finite Element Method (FEM). For this, I am assuming that a scalar plane wave propagating in the z-direction is illuminating a 1D screen with an aperture at the center. To set up the model, I am assuming the following Neumann conditions for the dependent variable u(x,z):

Absorbing boundary condition (ABC) to model an open system (boundaries at infinity) and ABC for the diffraction obstacle too, $$n\cdot \nabla u=-ik_0 u$$

radiation condition for the aperture to model illumination by a plane wave, $$n \cdot \nabla u= ik_0 ((1+n\cdot e)-u)$$

where $$n$$ is the boundary normal and $$e$$ is the incoming wave propagation direction. The solution obtained with this model is similar to the one using the KDI, but there is a mismatch in the scale and phase of the solution

I have tried different boundary conditions, for example, Sommerfeld's and Kirchhoff's, the ones used to deduce the KDI, but the results are similar. The solutions differ in a factor of $$\sqrt{s}$$ in the 2D case.

I found that the numerical and theoretical solutions match in the 3D case, i.e. u(x,y,z) and the aperture being a surface instead of a line as in the 2D case. The problem remains in how to reduce the KDI to 2D instead of 3D, normally when one consider a 2D diffraction problem, one assumes that, for example, the aperture extend to infinity in the $$y$$ direction, therefore it is not necessary to take the $$y$$ coordinate into account. I think that maybe the $$\sqrt{s}$$ factor comes from the reduction to 2D of the KDI. I tried doing the integration from $$[-\infty,\infty]$$ numerically, but it didn't work.

I will appreciate if someone knows what is wrong with the model or how to reduce de KDI to 2D rigorously.

Note: The boundary conditions were taken from Jian Ming's "The Finite Element Method in Electromagnetism".

• A guess: in 3D, the electric potential goes as -1/s; but in 2D, it should go as Log(s). Should come from Gauss divergence theorem/Helmoltz eq , hence the Log(s) should appear in the Green's function. Put the Log(s) in the KDI instead of 1/s? Needs some trickery as elecromagnetism in 2D is divergent; your attempt to do the integral in y probably crashed against this divergence. – patta Apr 12 at 6:49
• You may try to integrate the KDI / (2y) from -y to +y; then let y go to infinity...? – patta Apr 12 at 6:58

Short instruction. In 2D case one needs instead of $$\frac{e^{ikr}}{r}$$ to put a Bessel function of the second kind. The 2D version of KDI is the following: $$u(r)=\frac{1}{4} \left[ \int_S u \partial_n Y - Y \partial_n u \right]\quad(1)$$ where $$Y=-Y_0(k r),$$ where $$Y_0(x)$$ is a Bessel function of the second kind $$\pi Y_0(x)=\int_0^\pi \sin(x \sin{\theta})d \theta - \int_0^\infty e^{-x \sinh(t)}dt.$$ For large $$x$$ one can use an approximation $$Y_0(x)\approx-i\frac{1}{\sqrt{2\pi x}}e^{i(x-\frac{\pi}{4})}.$$ You can notice that your phase shift is about $$\pi/4$$ everywhere, so this formula is quite likely applicable to your case.
Some background. For any two functions, $$u$$ and $$u_1$$, satisfying the Helmholtz equation, the following integral can be shown to be zero $$\int_S u \partial_n u_1 - u_1 \partial_n u = 0.\quad (2)$$ Function $$u_1=Y_0(k r)$$ satisfies Helmholtz equation everywhere except at $$r=0$$. Picking an area that includes $$r=0$$, and cutting out a tiny disk out around $$r=0$$ ensures that the integral is zero over such area. The integral over the perimeter of the tiny disk gives us $$u(r=0)$$ (with some coefficient), while the integral over the big outside perimeter gives us the right side of the Eq. (1).