Feature request: gyrotropic media

[FilipDominec]

I would like to re-open the topic of gyrotropic media in MEEP. (Related terms are also Faraday effect, nonreciprocal behaviour, magnetooptics etc. Mathematically, all these are represented by nonzero imaginary part of the off-diagonal components in the permittivity tensor.)

While I found no access to this functionality in its current version, there are hints it has been implemented years ago:

• "Complex dielectric tensors are also useful for gyrotropic (magneto-optic) media. We have some internal modified versions of Meep that support this. However, merging these modifications into an official Meep hasn't happened yet, partially because I have another student who has been rewriting the PML and timestepping code" (January 2008, https://www.mail-archive.com/meep-discuss@ab-initio.mit.edu/msg01383.html)
• "We are also grateful to (...) Y. Chong for early support in Meep for gyrotropic media." (http://ab-initio.mit.edu/~oskooi/papers/Oskooi10.pdf)

This would be a very interesting feature for many researchers. (I know it from personal collaboration with other people, and it has also been asked for repeatedly in the mailing lists.)

I believe the transfer of MEEP development to Github can finally get together people with the know-how with people having time and motivation to make it work: Could we merge the obviously existing nonreciprocal MEEP code, perhaps as a mere prototype or experimental feature?

[oskooi]

Support for gyrotropic materials would be a nice feature to add. This has been on the TODO list for quite some now. The work that was referenced previously was incomplete and never merged into the codebase. Currently, Meep only supports anisotropic, real-symmetric, permittivity tensors. To generalize this to Hermitian tensors (i.e., no absorption loss) for frequency-independent materials requires modifying the time-stepping routines, particularly update_eh.cpp, which involves computing E from D using ε^-1. We may also consider generalizing the susceptibility class to support frequency-dependent media.

[FilipDominec]

I would personally prefer keeping the core computation as it is, and introduce the purely gyrotropic perturbation as another callback - if possible. This way we will make sure that MEEP behaves exactly the same in terms of efficiency and numerical stability for most users, but the access to gyrotropic behaviour will be still easy and similar to the definition of Lorentzians.

[stevengj]

One question is what physics to implement here:

• The simplest thing would be a complex-Hermitian offdiagonal element in ε or μ, but this isn't really physical except for narrowband simulations, because it makes Maxwell's equations complex.
• The "right" think would be a nonsymmetric offdiagonal element in the polarization equation, maybe?

Some literature search would be appropriate to see what people have implemented for this in FDTD and what timestepping equations they used.

[stevengj]

cc @seewhydee, who has implemented this kind of thing in the past, in case he has any suggestions.

[seewhydee]

If memory serves, the code I implemented was based on method 1 in Steven's above comment: i.e., using the complex fields representation and directly doing the complex algebra during the D step. The code worked, but I have no idea where it is now.

For actual inclusion in Meep, probably the best way to implement this functionality is method 2: add a dynamic off-diagonal contribution to the polarization, which works with real fields. One downside is that some researchers just want a simple model of magneto-optic response without worrying about dispersion. But it has the virtue of consistency with Meep's existing handling of anisotropic media.

Not sure about the details of the actual timestepping equations. Will dig around in my notes and let you know if I find anything.

[seewhydee]

Looking a bit closer, I think a straightforward approach is to tweak Meep's existing polarization model for Drude-Lorentz susceptibility, by adding an extra term γ dP/dt × n to the left-hand side of the equation of motion for P. Here, γ and n would be a user-defined cyclotron frequency and unit bias vector respectively.

It should be possible, with minimal alterations, to add this term to the update_P routine in susceptibility.cpp. It's an energy-conserving precession in P, and thus seems unlikely to cause instabilities or any other funky issues. At a steady state frequency ω, this term would yield a contribution to the susceptibility tensor with Voigt parameter γ / ω.

As mentioned, this approach might be annoying to people who just want to stick specific real on-diagonal and imaginary off-diagonal values into epsilon. They'd have to figure out the right combination of cyclotron frequency, operating frequency, and even dispersion parameters like ω_n. And since dispersion is intrinsic to the model, the desired epsilon tensor would only work at a single operating frequency. But... it would be consistent with Meep's existing handling of complex materials.

[stevengj]

The complication in implementing this is that, right now in update_P the d²Pᵢ/dt² equations are completely decoupled for each component of Pᵢ. Adding an off-diagonal dP/dt term requires this to be rewritten — especially since dP/dt involves both the new polarization and the old polarization, so you have to solve for the update equations differently.

Yidong, do you have a reference on what the revised finite-difference update equations look like?

[seewhydee]

meep_gyrotropy_note.pdf

I don't have an independent reference handy, but here are some notes. It seems that the update can be done if the p -> pcur -> pp overwriting isn't immediately performed for each component. In other words, the polarization at time step t+1 could be saved to a pnew slot, so that p and pp are still available for the update calculations of the other components. Then after all components are done, p and pp can be overwritten.

Or it could be done using some kind of split-step updating, I dunno.

[stevengj]

Not immediately doing the overwriting seems like it would require another temporary array.

Possibly the best option would be to simply update 3 components at once in the loop?

[seewhydee]

The main problem I see with updating 3 components simultaneously is that LOOP_OVER_VOLUME_OWNED has a component (c) input. So the loop over components/directions needs to be outside LOOP_OVER_VOLUME_OWNED, instead of the other way round.

Using a temporary array also makes it possible to stick quite closely to the existing Lorentzian susceptibility update code.

[seewhydee]

Here's a trial implementation. C++ code changes in gyrotropy-1.patch, Python/Scheme stuff in gyrotropy-2.patch. Not properly tested yet, comments on the overall approach/strategy welcome. (I may be able to get an undergrad student to help with testing.)

gyrotropy-patch.zip

[seewhydee]

And here is my working for the matrix inversion code blob: meep_gyrotropy_note.pdf

[stevengj]

Can you make a github pull request rather than posting a zip file? It is much easier to review patches in a PR.

(You can also eventually put documentation in an .md file with the PR in the docs/ directory — these can include LaTeX equations $...$ and $$...$$.)

[stevengj]

Yeah, I guess it's hard to avoid temporary arrays due to the averaging required for off-diagonal components, similar to computing E from D.

Is there much need for a gyrotropy vector b that is not constant? There's always the question of what happens at material boundaries (σ goes to zero outside of the material) … not handling this properly can lead to instabilities.

[seewhydee]

OK, I've opened pull request #863.

I'm not 100% sure what you mean regarding "gyrotropy vector b that is not constant". If you're referring to having b varying continuously in space across a given material, sure, that would be nice for modeling the effects of continuously varying magnetic fields. But based on my limited understanding of the code architecture, it seems like a lot of work to implement, just for the sake of a relatively uncommon scenario. In a typical setup for (say) topological photonics, each domain would have a uniform gyrotropy vector, representing a ferrite magnetized in a set direction, though different domains can be magnetized in different directions.

(If I understand right, there's a similar situation with the multilevel atom code. It would be nice to have continuously varying atomic parameters, to model spatially gain media; but it's not a deal-breaker.)

Regarding material boundaries: we'll see if testing turns up any instability problems, but it would surprise me since this is such a minor small tweak to the polarization dynamics, and the polarizations are damped anyway.

One aspect of the current setup that does seem a bit ugly is that the different susceptibility types are handled completely separately. So one can't specify a noisy gyrotropic medium, or a gyrotropic medium with multilevel atom dynamics, etc. Seems like a kinda arbitrary limitation.

[stevengj]

One aspect of the current setup that does seem a bit ugly is that the different susceptibility types are handled completely separately. So one can't specify a noisy gyrotropic medium, or a gyrotropic medium with multilevel atom dynamics, etc. Seems like a kinda arbitrary limitation.

It's a tradeoff. Having them separate simplifies the implementation and the testing, but on the other hand you're correct that it limits the kinds of physics that can be modeled.

Since noise, gyrotropy, and multilevel atoms are all kind of esoteric features, though, I'm not especially concerned that you can't have them simultaneously in the same polarization state right now, though — it's worth it for simplicity/orthogonality of the implementations.

[mawc2019]

Is it valid to set bias in mp.GyrotropicLorentzianSusceptibility() as an imagery vector, such as something like bias=(0,0,1j)? No error message is reported when I do so, but the results are not expected. @stevengj

[stevengj]

No, I don't think so.

[mawc2019]

My purpose was to obtain a real but asymmetric permittivity tensor, so that the system is nonreciprocal but has time-reversal symmetry. As bias has to be a real vector, the other way to get the desired permittivity is to set sigma as an imaginary number, frequency as the source frequency, and gamma as a nonzero real number. However, simulation fields are NaN or Inf is reported. Perhaps mp.GyrotropicLorentzianSusceptibility() could not be used to generate a real but asymmetric permittivity tensor?