Under which circumstances is Snell's law applicable?

Question

For a few years I have been tinkering with a numerical solver for light propagation. I stuck with the Helmholtz equation in the paraxial approximation because of lasers.

$$\partial_z f= \frac{i}{2k}\nabla^2 f.$$

Lately I have been trying to use that PDE to model light propagation through a medium that has a variable refractive index with respect to the spatial coordinates. I have started from the Helmholtz equation

$$\nabla^2 E + k^2 E = 0$$

and then applied the same ansatz for the solution

$$E=f\exp(ikz)$$

with the main difference that now $k=k(x,y,z)$ since $k=2\pi n/\lambda$ and $n$ is a function of space.

I have implemented a numerical solver for computing $f$ from the resulting PDE and I tested it with what I think is the most basic scenario I could think of, namely refraction. So I defined the refractive index by splitting the spatial domain in each region and I did the computation. For numerical reasons, I had to apply a smoothing function where the two regions meet in order to keep my numerical approximation from blowing up.

What I found is that refraction does indeed occur, but the angle varies with the width of the blurred area between the two regions. Testing this with an online geometric optics simulator I got kind of the same result, namely that by creating thin layers with an increasing refractive index, the angle at which refraction occurs changes.

This got me looking for explanations or cases in which others might have tried what I did and I did not manage to find much. Either you get the refractive index directly from manipulating the Maxwell equations (as I found in Principles of Optics by Wolf and Born), or some handwavy links people do between Helmholtz, plane waves and conditions that can be imposed at the region where the two media meet, or the eikonal equation which gets into ray optics territory where Snell's law is applied.

So, is Snell's law a result that is good enough for most practical purposes but it does not stand as well if one goes into the wavy nature of light or is my approach just wrong because I introduced the varying refractive index somewhere in between my equation and Maxwell?

EDIT:

This is a plot of a beam that propagates from the upper border to the bottom one, passing through a medium that looks like the one in the second image. In the refractive index profile, purple is 1 and yellow is 1.33.

The steps on each axis are $10^{-8}$m, I know, not physical, but it is for the numerical solver. Since the dimensions are the same on both axes, I computed the angles and it was close. I did the computation again for the angles and I got $i-r = 0.1026 rad$ from the solver, and $i-r = 0.1207 rad$ from the geometry and Snell law. It's close, but only because I chose a smooth parameter that got me a value that is close to Snell.

Answer 1

Snell says that in a layered medium $n(x) \sin\theta(x)= {\rm const.}$, so, in the domain of validity of ray optics, the angle of incidence $\theta(x)$ should only depend on the current value of $n(x)$ and not on the paticular profile that interpolates between two different media. It is hard therefore to see how your geometrics optics simulator can disagree with this.

In the wave optics case, the refraction angle should not depend the vector character of Maxwell's equations: a scalar field should give the same ray path --- although the refelection coefficients will differ. The paraxial approximation does not handle reflection though. Do you have some some plots of your output?