update on multilevel atoms

[acerjan]

So I've been working through Ardavan's multilevel atom implementation, I've found at least one bug (listed below), and I think I'm getting close, though I still haven't gotten laser thresholds to agree correctly. (I expect to be able to work through this, but i'll update if I stall out.)

However, I did want to solicit a bit of help with one change -- in lines ~250-310 of src/multilevel-atom.cpp, the code is updating N(t), the populations of the atomic levels, and this is dependent upon the electric field and polarization. In the conventional writing of the oscillator model of the gain atoms, this coupling takes the form of E dP/dt. However, for perfect agreement with SALT, there is an additional term which is proportional to E P. (this was part of the punchline from the notes I sent out a while ago.) To get this term, can I just copy lines 287-300, but replace instances of "dP = p - pp" with "Pave = p + pp", but divide by 2 and multiply by dt? I assume I'll still need to divide by 32 as well?

The previous bug(s) were (a) in the definition of alpha in scheme/structure.cpp, around line 1280. This is 1/(hbar omega), with hbar = 1. However, as frequencies were being specified without the 2pi, alpha needed to have the 2pi added back in by hand; and (b) there was some funny business going on with sigma in line 361 of multilevel-atom.cpp.

[oskooi]

Here is a reference for others to follow this discussion involving the equations implemented in #273.

[oskooi]

σ in the polarization update equation (number 2 from the notes figure above) is actually a representation of two different sigma terms in the code: one for the overall susceptibility class and the other for the particular transition of the multilevel system. That's why in multilevel-atom.cpp#L356 σ appears as st * s[i]: st is a meep::realnum *sigmat = new meep::realnum[T * 5]; defined in structure.cpp#L1314 for each two-level transition whereas s[i] refers to realnum *sigma[NUM_FIELD_COMPONENTS][5] of the susceptibility class. Unfortunately, the notes did not make this clear.

[acerjan]

Alright, I've achieved validation. Here is an updated plot comparing MEEP to my existing FDTDs and SALT for the one-sided n=1.5 Fabry Perot cavity. And here is another comparison very close to the laser threshold. As you can see, the agreement between MEEP and my FDTD using the "uncorrected" oscillator model equations is quite good. For the 3rd lasing mode, in green, I suspect I didn't run MEEP quite long enough to resolve the low intensities perfectly. Modes close to threshold require much longer runtimes to properly resolve, so the mild discrepancy here doesn't bother me.

Ok, so this means all we need to do for agreement with SALT is to "fix" the oscillator equations in MEEP, just as we discussed before, by adding back in the two terms usually dropped in the derivation of the oscillator equations from the bloch equations.

[stevengj]

Thanks again for working on this!

[stevengj]

I should mention that Song Liang Chua deserves much of the credit for the initial draft implementation, which he developed sitting together with me back in 2010. Song left and discontinued work on it, so the code was subsequently sent to Ardavan, who checked it into the repo and began debugging and further development.

[oskooi]

Congrats Alex. It's satisfying to know that all the hours we spent discussing and reviewing the derivation and implementation for this feature during the last two years was not in vain. Looking forward to merging your updates, adding a Python frontend, and using it to investigate lasing in complex structures.

[acerjan]

Oh no, of course, Song and Ardavan deserve most of the credit for this. Really, this was all in excellent shape except that the code's alpha needed a factor of 2*pi, and sigma was being inserted as sqrt(sigma), which wouldn't stop the system from lasing. The rest of the validation was simply carefully going through all of the SALT units, which are unnecessarily confusing. All of this is now laid out in the heavily revised example code, multilevel-atom.ctl.

Finally, I've altered the algorithm MEEP is using slightly, to add back in the two terms which are usually approximated to zero in the oscillator equations, and now the results nearly perfectly agree with SALT, and my previous Bloch FDTD. And close to threshold: There might be some slight funny business happening at D0 > 0.6, the code appears to either have a 4th lasing mode on or some light 3 wave mixing which is dumping some intensity there that would go away with longer runtimes, and I'll check this later today, but I do think we're in good shape now. I'll set up a pull request.

[oskooi]

As a useful reference for this feature in addition to the tutorial example (which we need to port to Python once the frontend is ready), we should also create a new "Saturable Absorption" page under the "Features" tab of the documentation. This will be used to provide a brief theoretical background of modeling multilevel gain media in Meep including a discussion of the actual oscillator equations which are implemented in the code. We should also include the figures posted here which demonstrate agreement with SALT and your Bloch FDTD solver. Note that the documentation is part of the repo as markdown files in the doc/docs.

[acerjan]

Ok, i can take a look at this in a couple of days.