Implementation of gyrotropic susceptibilities #863
Conversation
|
Overall, this looks great! Have you tried it out on any example problems yet? Once there are some minimal tests that indicate that this is essentially working, I would be happy to merge as an experimental feature since it doesn't really affect the rest of the code. If you have an example reproducing something from a paper to use a tutorial that would be good too. My biggest concern (modulo bugs) is whether this will introduce instabilities, especially at the boundaries of gyrotropic media since inhomogeneous anisotropy is always a bit tricky in a Yee grid. |
|
Haven't done any rigorous testing yet, will work on that soon. Off the top of my head, there are a few straightforward tests to be done in 2D:
For 3D, I'm not sure; it may be necessary to compare to various Comsol simulations. |
…ady in gyro_tensor.
…ghly more stable...
|
Sorry for the long silence. Here's an update. I wasn't too happy with the equation of motion that just adds a precession term to the Lorentzian damped harmonic oscillator, since it does not reproduce the right frequency dependence for materials like ferrites. So I spent some time implementing another equation of motion, the Landau Lifshitz Gilbert equation, which does give the right frequency dependence (this is on the "landau-lifshitz-gilbert" branch). Implementation notes here. However, there's a stability problem that I haven't been able to fix. For relatively weak gyrotropy and short runtime, the simulation behaves as expected. Here is a figure plotting Ex and Ey for a wave propagating along a gyrotropy axis z:
As desired, the Ex component precesses into Ey. But with increasing t, a spatial instability sets in and creates an oscillation on the length scale of the grid. It doesn't seem to be a problem with the Euler integration of the Landau Lifshitz Gilbert equation by itself -- for an isolated point, the integration doesn't seem to exhibit any instability. There seems to be some kind of bad interaction with the EM fields on the grid. Kinda stumped right now, so suggestions welcome. |
|
I also found a bug in my code (forgot to change the point source to be a line source when I added a finite thickness), but even with the newly updated version, my bandgap is still a bit smaller than I expected. I've also been through your new set of notes and I agree with your derivation. It may be that the smaller bandgap I'm seeing is entirely correct, and is due to changes in the effective dielectric tensor for frequencies further away from the central frequency, and that the thing to do is run a series of very narrow sources to build up the band structure that way. For reference, the bandgap I'm expecting is ~20% of the central frequency. This change will take a bit of time to implement, so it'll have to wait until I return from holiday. |
|
@acerjan: the topological band gap you're testing on is really wide, and dispersion is highly likely to be a significant issue. For the quantitative verification, I suggest simply testing against a ~2-3% bandgap. @stevengj: I think the code is converging, so could you take a look at the proposed interface and comment on whether it's acceptable? (The Scheme interface will track the Python one, once finalized). |
|
@seewhydee, the interface looks good. |
|
If we can have a sanity test of this feature, I think we can merge it? |
|
Fine by me. As of now, Faraday rotation of a plane wave quantitatively matches theory for all three types of gyrotropic media. For inhomogenous structures containing both gyrotropic and non-gyrotropic media, the simulations run without field instabilities, but we have not quantitatively checked the results yet. As a clearly-labelled experimental feature, it's OK to merge whenever you like. Alex will compare the topological photonic crystal band gap size when he gets back, and there are a bunch of other tests I'm interested in doing (e.g. looking at the dispersion relation of magnetic surface plasmons) before removing the "experimental" label. |
Great! Can you add an automated test to |
|
Looks good to merge once tests are green. |
|
The |
|
Oh, right. Added it. |
|
After returning from a couple of different trips, I've finally had a moment to run simulations which account for the dispersion in MEEP's model of gyrotropic material relative to MPB's. Although Steven's suggestion of doing self-consistent frequency calculations in mpb would probably have been less computer runtime, just running a bunch of meep simulations with different choices of the gyrotropic media was faster to set up, and computer resources are bountiful. As a reminder, I'm simulating the boundary between a topological and a trivial photonic crystal, so we expect to see a chiral edge state. Of course, mpb finds both edge states, but meep should only find one, as the other border is buried in the absorbing region. As can be seen in this plot comparing mpb (red lines) with meep (blue +s), the results are similar, though some slight differences are seen -- the meep frequencies seem to be reduced by ~2-3% compared with mpb, but the overall shape, and bandgap are nearly the same. Is this a typical amount of discrepancy between the two simulation techniques? Or is this indicative of an error?
If you think the error is small enough to not be concerned, I'm happy to prepare some of this as an example for meep / mpb, just let me know what you'd like to have added. |
|
Could this be a due to discretization? Does changing the resolution alter the frequency shift? It certainly would be great to have a photonic TI as a Meep tutorial, but doing repeated simulations over a frequency sweep seems like an unnecessary complication for that purpose. Could you try picking parameters where the bandgap is about 1/10 the present size (say a 3% bandgap rather than a ~20% bandgap)? Then dispersion effects ought to be much less pronounced, and it should be possible to compare Meep results to MPB results (using mid-gap values) without a frequency sweep. |
|
Yeah, i'll move to a 3% bandgap or something smaller so i can check the discretization more easily. It still might require 2 or 3 runs, but not the ~26 used in the above plot. |
* Implement gyrotropic susceptibility class. * Add Python and Scheme support for gyrotropic media. * Initialize bias vector in python susceptibility struct. * Remove "bias" from gyrotropic_susceptibility; the information is already in gyro_tensor. * Add gyrotropic media to docs * Minor copyedit * Fix logic in py_susceptibility_to_susceptibility * In add_susceptibilities, always pass a 3-vector as gyrotropic bias * First try at gyrotropic media tutorial * Fix errors in gyrotropy formulas in notes * Re-implemention of gyrotropy using LLG equation * Tweak handling of sigma tensor in gyrotropic case * Update Python doc. * Drop 2pi factor from gyrotropic sigma * Doc updates * Fix last change to update_P * Fix printf typo * Remove spurious 2pi factor in alpha (which is not a rate) * Minor code tweak * Use a central-difference scheme for the LLG dynamics, which seems slighly more stable... * Try implementing the full nonlinear LLG equation * Add implicit static polarization to gyrotropy implementation * Put static P back in in subtract_P * Add gyrotropy example * Fix; use LOOP_OVER_VOL instead of LOOP_OVER_VOL_OWNED to ensure updating of off-diagonal components * Clamp the magnitude of the LLG polarization vector. * Revert inadvertent unrelated change to meep.i * Minor code cleanup * Flag "needs_W_notowned" for gyrotropic media * Update gyrotropic P components explicitly; don't use LOOP_OVER_VOL_OWNED * Enable needs_P on all components for gyrotropic media * Fix gyrotropy scheme to track 9 polarization components per unit cell. * Revert unrelated last change to meep.i * Avoiding need for allocation of P_tmp in gyrotropy_data. * Implement num_cinternal_notowned stuff for gyrotropic media * Update documentation for gyrotropic media, and relax some minor restrictions. * Add virtual keywords to gyrotropic_susceptibility methods * Merge latest changes from master * Remove gyrotropic-dispersion.py (incomplete attempt) * Complete merge * Update Materials.md to discuss both Lorentzian and LLG gyrotropic models * Introduce a new gyrotropy_model enum type, to allow for the LLG model. * More plumbing to provide support for Landau-Lifshitz-Gilbert type gyrotropy model * Fix typo in susceptibility update equation * Fix typos in Faraday rotation formula in docs * Merge from master * Reimplement linearized-LLG updating equations * Fix typo * Minor code clarification * Fix Faraday rotation example * For Landau-Lifshitz-Gilbert model, ignore the magnitude of the bias vector. * Fix minor hiccup in docs. * Support dumping and undumping of gyrotropic susceptibilities * Doc updates and minor tweaks accompanying last merge * Translate Faraday rotation tutorial from Python to Scheme * Fix typo in last change * Fix missing 2pi factor in gyrotropic LLG susceptibility's sigma parameter * Minor fixes for gyrotropy documentation * Add Faraday rotation unit test * Use absolute tolerance (in degrees) for Faraday rotation unit test * Add faraday rotation test to python/Makefile.am
Closes #60 ...