FDTD Modelling of dipole emission in 2D materials


Sorry to bother you. I had a couple of questions about modelling dipole emission in 2D materials using lumerical. I have been playing around a little with some simple dipole simulations - but I have been struggling a little to understand how to set things up so that I am convinced that lumerical is even noticing the very thin layer of 2D material - (ie so the mesh is small enough around it) - without making the simulation time so long as to make simulating the more complicated structures I actually want to model difficult.

I guess that the actual material properties of the TMDCs might cause a problem as I go down to small sizes - but for the moment I’m happy with fudging the material database.

I am currently attempting to model a hexagonal boron nitride layer on various materials eg SiO2 on Si- and trying to calculate the collection efficiency - and the quantum efficiency/Purcell factor. Since I’m aiming simulate structures later which will have Purcell enhancement - and since I’m going to be overiding the mesh size to force lumerical to see the hBN I guess I may need to use transmission boxes to calculate the dipole output (rather than dipolepower as discussed in various other posts) - Purcell factor: dipolepower vs. transmission box

I have arbitrarily settled on a hBN thickness of 4nm (which I think (?) in terms of the material is roughly 10 layers) - though probably being able to vary the thickness would be a good idea.

My question(s) at the moment are:
If I want to create a transmission box - does it have to be all inside the hBN? - If so that would mean that the mesh size should be at most 0.2 nm to give 10 mesh boxes inside the Tbox. I guess if I can allow the T box to extend outside the HBN then it would allow me to use a coarser mesh.
Does the T box need to be a cube? ie. Can I get away with a larger TBox in X and Y than in Z - and then also be allowed to have larger mesh overide size in the X/Y directions.

I noticed that there were some layered stack analytical tools which might be helpful as well for this project - perhaps doing FDTD on such thin layers isn’t ultimately a very sensible way of modelling things - though I think they are not on the standard licence. I was wondering whether it was possible to arrange to test them out (and maybe get an idea how much extra they are) - I think our current set of licences are up for renewal in about a month.

Hello @HaddenJ,

Thanks for posting on this fascinating topic and welcome to the community. While I do not have experience modelling 2D hBN it seems that it is an insulator with optical properties weakly dependent on it’s thickness; therefore, I would not expect it to present too great a challenge.

As long as it is much thinner than the wavelength of light it should appear 2D to the radiation. I would use a mesh override to force a single mesh cell to completely cover the hBN layer, and experiment with varying the thickness to see if it changes the reflection and transmission properties greatly. That is change the mesh override extent and thickness at the same time. I do not know the answer, so convergence testing is your best bet, but I would suggest starting larger, $ \lambda /10 $ for example, and decrease that by a factor of 10 to see if your results change.

I am not sure what sort of published values may exist for the resistivity, but you could possibly use a 2D conductive (non-dispersive) or 2D sampled data material (dispersive). These materials are implemented differently, but I would still recommend using a mesh override region.

It would be a good idea to use a transmission box as a sanity check, but it should only really differ from the purcell factor returned by the dipole if their is loss in the material. The transmission box does not be contained within the hBN, and can be any rectangular dimension yet monitors need to align with the FDTD mesh for accuracy.

The STACK solver is a frequency domain technique and can be very useful for large stacks with high Q. I think it might be helpful for your proposed workflow, for comparison and validation.

In regards to licensing I would ask that you get in contact with sales@lumerical.com, they may will be happy to assist you with any licensing questions and requests you may have.


Hi there @trobertson,

Many thanks for your reply. Could I just clarify - so I’m understanding you correctly:

As far as setting up the geometry - you are saying essentially I can increase the grid size of my mesh so it roughly matches the thickness of my HBN layer?

I will need to tweak the mesh overide a little so that the dipole which sits within the hBN is as usual halfway between two nodes of the mesh - so for a x dipole it will likely need to be 2 mesh boxes thick, and for a z dipole it will need to be 3 mesh boxes thick.

I can use a transmission box as a sanity check - but it doesn’t need to be completely inside the hBN.

Regarding the material properties - I guess I will need to play around a bit - for the first scratch attempt I was treating it as a dielectric with the experimental refractive index. (simply because I wasn’t sure what to do - and I’ve previously mostly modelled dielectrics).

Hello @HaddenJ,

It would be best to make it so that an integer number of mesh cells cover the hBN and the mesh cell boundaries align with the interfaces. That means that the mesh override region exactly matches the thickness of your 2D material. Then test how thick you can make the region by starting with a large value and shrinking it until you stop observing changes i.e. convergence testing.

I am not sure what the behavior will be when trying to place a dipole source within a “2D” material. It will likely seem to be a cavity if it is given finite thickness, and so you would see some Purcell enhancement. Whether this is physically realistic I do not know. It will not be possible to place a source inside a 2D material should you choose to define it that way.


Hello @HaddenJ,

I am not sure if you have seen this post, but it involved modelling a 2D material

The very small override region in 3D which we discussed should give you the correct answer, but it may cause some other issues with the PML. I now believe that the 2D material model may be a better approach.

You can use the script command stackrt amd stackfields without a license to test your material models.


Thanks for your help. I have been playing around a bit more with my FDTD and stackrt simulations - and reading through the examples you sent me.

I have a quick question about the stackfields function and how to understand the ‘fields’ which are returned-when theta is not zero. Apologies in advance if this is a silly question - I feel like what I’m guessing here is correct - but it isn’t exactly clear from the stackfield description what it is actually returning - and I’d like to make sure I understand it…

ie the array of 6 dimensions which describes Es and Ep. Where I’m talking about Ex,Ey,Ez here I mean the 6th dimension when you return result.Es and result.Ep.

say you calculate:
result= stackfield(index,d,f,theta);
When θ = 0 - the trivial case - everything seems clear enough
Es appears to have all it’s field in Ey - so that it makes sense to use Es = pinch(result.Es,6,2)
Ep appears to have all it’s field in Ex - so that it makes sense to use Ep = pinch(result.Ep,6,1)

When θ is not zero:

it seems like the field for Es remains the same - so that Ex and Ez are both = 0,
but for Ep, the field is shared between Ex and Ez. I guess that in this case Ex=Ecos(θout); and Ez=Esin(θout); - where θout is whatever the snells’ law tells us the plane wave’s theta should be throughout the stack.

I guess that I had been wondering whether it is possible to keep track of the angle of the light going through the stack - but presumably it is a case of calculating it from the Ex,Ez at each point so theta=atan(Ez/Ex)

Is that correct?

Finally - the result suggests that the Es, Ep values are returned vs x,y as well as z, freq/lambda, and theta. I assume that the x,y are essentially always redundant - since it is a 2D stack - or is there somethign I’m missing?



Hello @HaddenJ,

The P polarization refers to light in the plane of incidence, defined as containing the k vector and the normal vector to the interface. Thus with STACK this plane given as the xz plane. S polarization being normal to the plane of incidence will always be in Y. I had to double check the results, but I believe you are correct, for all the relations that you suggest.

To keep track of the k vector of the plane waves it would be straightforward to perform a loop over the index vector, calculating $ \theta_t $ via snells law for each layer which would then become $ \theta_i $ in the next.

Since these are plane waves, the xy coordinates are redundant since it is assumed to be uniform in this plane. I suppose these were kept avoid ambiguity and make it easier to concatenate multiple matrices later.

I hope this helps.