Is my simulation result for unpolarized light correct?

This is a follow-up of this question. After that, I picked up some knowledge of FDTD (an algorithm for solving Maxwell's equations) and simulated following scene:

Pic 1

As the picture shows, a silicon cone ($r=5 μm, θ=30°$) is placed in vacuum (truncated with PML). Incident light ($λ=532 nm$) is perpendicular to the bottom of the cone. In my simulation the incident light is modeled as a electric field hard source placed at the bottom of the cone and its electric field intensity obeys Gaussian distribution i.e.

$$E\left(x,y,z_0\right)=A_0 \sin (2 \pi f t) e^{-\frac{x^2+y^2}{w^2}}$$

Here $E$ can be either $E_x$, $E_y$, $E_z$, $z=z_0$ is the plane at which the bottom surface of the cone is.

The whole domain is discretized with 200×200×400 uniform grids ($\Delta x=\Delta y=\Delta z=50 nm$), along with 20-grid-thick PML layer enclosing it.

Here's my simulation result of the light intensity ($|E|^2$) at the plane $100 nm$ away from the tip of the cone for z-polarized light (notice that the z-axis is parallel to the central axis of the cone in my coordinates):

Pic 2

and y-polarized light:

Pic 3

So here comes the confusion. As Steve B suggest in his answer, "the intensity for unpolarized light is the average of the intensity in the two simulations (for polarized light)", but since this two result is… not independent of the rotational orientation of polarization analyzer (if there exists one), a simple average of them won't be, too:

Pic 4

This is the average of z-polarized and y-polarized result, it's almost the same as that of the y-polarized light because that of z-polarized one is in fact very weak.

So the unpolarized light is no longer unpolarized after it passes the cone?

Moreover, it's not hard to realize that simulation result for x-polarized light is just the same as the y-polarized one, only a 90° rotation needed, so the average of x-polarized and y-polarized light is:

Pic 5

Still, the picture don't have full rotational symmetry, what's more, it's different from the average of the z-polarized and y-polarized one! So I'll get different simulation result for unpolarized light with different average method?

Is my way to "take the average" wrong? Or it's just the truth?

My code isn't contained in this post for its lengthiness. I think it's trivial since it has been checked for times and nothing seems to be wrong.

Any help would be appreciated.

Wrong Guess Update 1:

In the comments below, DumpsterDoofus and Ruslan guessed that the PML I used in the model maybe not work as expected i.e. the lack of full rotational symmetry of my image may be caused by the reflection of the boundary I use.

If this guess is correct, the image I get will be different when I change the size of simulation domain, so I tried two more simulation, and the result, in a word, proves negative.

To make the simulation faster, I've chose a smaller cone ($r=\frac{5}{2} μm, h=\frac{5 \sqrt{3}}{4}μm$) this time. Two domains with different sizes (200×200×100 grids and 120×120×70 grids) are used in the two simulations separately. Another tiny change is that this time I've turned to additive source instead of the hard source used in the simulation above. Grid size and the thickness of the PML layer (for this part, see my edit above) remains the same as that of the simulation above.

Here's how the whole domain looks like in the two simulations (the sizes of the cones are same).

Domain formed with 200×200×100 grids:

Pic 6

Domain formed with 120×120×70 grids:

Pic 7

And here's the result (the picture is still taken at the plane $100 nm$ away from the tip):

y-polarized light (simulation result of the larger domain):

Pic 8

y-polarized light (simulation result of the smaller domain):

Pic 9

Supposed unpolarized light (simulation result of the larger domain):

Pic 10

Supposed unpolarized light (simulation result of the smaller domain):

Pic 11

Apparently there is no significant difference except the spot size(Note: I forgot to adjust the DataRange of ListDensityPlot 囧, and I decide not to modify it in the following part of this post), and we still don't have full rotational symmetry.

Wrong Guess Update 2:

Ruslan then suggested that the round source I chose may generate unequally widening wave in x and y direction depending on polarization i.e. it maybe not work as expected.

If this guess is correct, the simulation result won't show perfect symmetry even if I remove the cone i.e do the simulation in an empty space. So I simulated the y-polarized light in vacuum without the cone, and the result shows perfect symmetry, so the source seems to work as expected.

The following picture shows how the whole domain looks like in this simulation (formed with 120×120×70 grids):

Pic 12

Here's the result (taken at same place as pic 9):

Pic 13

-
What would first surprise me is not lack of any rotational invariance of the intensities, but the difference of pictures of intensities. If I understand you correctly, the difference between two images should just be rotation of $\vec E$ by $\frac\pi2$, because you just rotated the source around cone axis, and the medium has symmetry with respect to rotation by arbitrary angle around this axis. But what you present is completely incompatible with this. How does your picture change if you continuously rotate the polarization angle from $0$ to $\frac\pi2$? – Ruslan Apr 1 '14 at 11:30
@Ruslan The difference between the two image is indeed a rotation of $\overset{\rightharpoonup }{E}$ by $\frac{\pi }{2}$, but it's not rotated around cone axis, z-axis in my coordinates is parallel to the cone axis while y-axis is perpendicular to it. My original description seems to be a little obscure, edited :) – xzczd Apr 1 '14 at 11:52
Hmm. How could light in vacuum be polarized in the direction of propagation? – Ruslan Apr 1 '14 at 12:01
@Ruslan !(Check the wikipedia…) Sorry for my poor basis of physics. Well, my original description isn't quite precise, in fact when simulating the incident light is just modeled as a electric source placed at the bottom face of the cone. By the way, the intensity for z-polarized light is in fact very weak, when compared with that of y-polarized light, it's almost masked. Still, even if I choose x-polarized and y-polarized light for my modeling of unpolarized light, my question mentioned above still exists. – xzczd Apr 1 '14 at 12:31
At this point your question has gotten so long as to be daunting. In my browser it spans 5-6 pages. Could you combine some of your figures into a grid and combine all of the updates into one coherent question? – Chris Mueller Apr 6 '14 at 14:17

As Ruslan said, your error lies in the fact that you used z-polarized light. There is no such thing as z-polarized light (it doesn't exist).

-
Sorry for my obscure description, see my comment above. – xzczd Apr 1 '14 at 12:37
@xzczd: You say that using the $x$ and $y$ polarizations, the problem still exists. I assume the $x$ polarized light looks like a $\pi/2$ rotation of the $y$-polarized picture you provided? I'm noticing that the pictures have the symmetry group of the boundaries (reflection symmetry). What happens when you widen the boundaries? It might be a boundary problem. – DumpsterDoofus Apr 1 '14 at 12:51
Yeah, your guess for the x-polarized light is right. The cone is modeled as being placed in infinite space, with the help of PML. – xzczd Apr 1 '14 at 13:19
Do you take PML as frequency-dependent? Gaussian wave, especially originating from a 2D source, is far from a monochromatic wave, so you may get reflection from a single-frequency PML. This may be the reason for lack of full rotational symmetry of your images. Try indeed to considerably expand your boundaries. – Ruslan Apr 1 '14 at 13:26
Gaussian beams have a decomposition into plane waves, and the $k$ vectors of those plane waves do not all point in the $z$ direction. Also, once the beam reflects off of a wall, there will be $z$ components to the $E$ field. – garyp Apr 1 '14 at 13:30

As mentioned by Ruslan , precisely speaking, what one should do to simulate the unpolarized light is to take an average of the intensity of all the orthogonal polarized light other than just 2 of them. Plane source is a special case because its z-polarized component is quite weak so it won't hurt even if only an average of x and y-polarized component is taken.

But wait, in Pic 4, OP has already taken an average of all 3 component, why he still can't get full rotational symmetry?

The answer is really simple but easy to ignore: the grid isn't dense enough, so it has formed a not-really-round cone.

After halving the grid size (with other condition remaining the same as Pic 9 and Pic 11), we got the following result:

y-polarized light (with same condition as Pic 9, except for a denser grid):

Pic 14

Supposed unpolarized light (with same condition as Pic 11, except for a denser grid):

Pic 15

z-polarized light:

Pic 16

Light intensity here is in fact very weak, the simulation result of passing rate is about 0.08%. BTW, grids for this simulation still seems to be not dense enough, but it's not a big deal.

To avoid any possible confusion, here is a more accurate simulation result of the z-polarized case:

z-polarized light (with same condition as Pic 16, except for a denser grid, $\Delta x=\Delta y=\Delta z=12.5 nm$):

Pic 17

Though still defective, the simulation result almost gains full rotational symmetry, being much better than Pic 16. The passing rate is about 0.04%.

-