Working and verified version of multilevel atom susceptibility. (#500)

* First commit by AWC to multilevel atom branch. Fixed 2*pi bug in the definition of alpha of /scheme/structure.cpp. Fixed a bug in how the E*(n2-n1) term in updating the multilevel polarization was working. Rewrote most of the multilevel-atom.ctl example to be more transparent with agreement with SALT. Still having trouble getting multilevel-atom.ctl to yield non-zero fields, I think this is a field initialization issue.

* Completely working version of the unmodified Oscillator Model Equations for gain media in MEEP. The results have been verified against my personal implementation of the oscillator equations in 1D.

* Working and verified version of multilevel atom susceptibility for lasing, including the two ~0 terms not normally written in the oscillator model equations. Candidate to be merged back into the master branch.

scheme/examples/multilevel-atom.ctl

@@ -0,0 +1,84 @@
+;; This file realizes a 1D, one-sided Fabry-Perot laser, as described in Fig. 2 of Cerjan et al.,
+;; Opt. Express 20, 474 (2012).
+
+;; Cavity definitions
+(set-param! resolution 400)
+(define-param ncav 1.5)                          ; cavity refractive index
+(define-param Lcav 1)                            ; cavity length
+(define-param dpad 1)                            ; padding thickness
+(define-param dpml 1)                            ; PML thickness
+(define-param sz (+ Lcav dpad dpml))
+(set! geometry-lattice (make lattice (size no-size no-size sz)))
+(set! dimensions 1)
+(set! pml-layers (list (make pml (thickness dpml) (side High))))
+
+;; For defining laser properties in MEEP, the transition rates / frequencies are specified in units of 2*pi*a/c.
+;; gamma-21 in MEEP is the Full-Width Half-Max, as opposed to gamma_perp, which is the HWHM in SALT.
+;; Additionally, the classical coupling element sigma = 2*theta^2*omega_a/hbar, where
+;; theta is the off-diagonal dipole matrix element.
+
+;; These different conventions can cause a bit of confusion when comparing against SALT, so here we perform
+;; this transformation explicitly.
+
+(define-param omega-a 40)                        ; omega_a in SALT
+(define freq-21 (/ omega-a (* 2 pi)))            ; emission frequency  (units of 2\pia/c)
+
+(define-param gamma-perp 4)                      ; HWHM in angular frequency, SALT
+(define gamma-21 (/ (* 2 gamma-perp) (* 2 pi)))  ; FWHM emission linewidth in sec^-1 (units of 2\pia/c)
+					         ; Note that 2*pi*gamma-21 = 2*gamma_perp in SALT.
+
+(define-param theta 1)                           ; theta, the off-diagonal dipole matrix element, in SALT
+(define sigma-21 (* 2 theta theta omega-a))      ; dipole coupling strength (hbar = 1)
+
+;; The gain medium in MEEP is allowed to have an arbitrary number of levels, and is not
+;; restricted to a two-level gain medium, as it simulates the populations of every individual
+;; atomic energy level.
+
+;; If you are using a 2 level gain model, you can compare against
+;; results which only simulate the atomic inversion, using the definitions
+;; gamma_parallel = pumping-rate + rate-21
+;; D_0 = (pumping-rate - rate-21)/(pumping-rate + rate-21) * N0
+
+;; In fact, even if you arn't using a 2 level gain model, you can compare against an effective
+;; two level model using the formula provided in Cerjan et al., Opt. Express 20, 474 (2012).
+
+;; Here, D_0 as written above is not yet in "SALT" units. To make this conversion,
+;; D_0 (SALT) = theta^2/(hbar*gamma_perp) * D_0 (as written above)
+
+;; Finally, note the lack of 4*pi in the above conversion that is written in many published SALT papers.
+;; This 4*pi comes from using Gaussian units, in which the displacement field, D = E + 4*pi*P, whereas
+;; in SI units, D = eps0*E + P, which is what MEEP uses.
+
+;; Gain medium pump and decay rates are specified in units of c/a.
+
+(define-param rate-21 0.005)                     ; non-radiative rate  (units of c/a)
+(define-param N0 28)                             ; initial population density of ground state
+(define-param Rp 0.0051)                         ; pumping rate of ground to excited state
+;; so for example, these parameters have D_0 (SALT) = 0.0693.
+
+;; Make the actual medium in MEEP:
+(define two-level (make medium (index ncav)
+			(E-susceptibilities (make multilevel-atom (sigma 1)
+			  (transitions (make transition (from-level 1) (to-level 2) (pumping-rate Rp)
+					     (frequency freq-21) (gamma gamma-21) (sigma sigma-21))
+				       (make transition (from-level 2) (to-level 1) (transition-rate rate-21)))
+			  (initial-populations N0)))))
+
+;; Specify the cavity geometry:
+(set! geometry (list (make block (center 0 0 (+ (* -0.5 sz) (* 0.5 Lcav)))
+			   (size infinity infinity Lcav) (material two-level))))
+
+;; Initialize the fields, has to be non-zero, doesn't really matter what.
+(init-fields)
+(meep-fields-initialize-field fields Ex
+			      (lambda (p) (if (= (vector3-z p) (+ (* -0.5 sz) (* 0.5 Lcav))) 1 0)))
+
+;; Specify the end time:
+(define-param endt 7000)
+;; Note that the total number of time steps run is endt*resolution*2. This is the origin of the extra
+;; factor of 2 in the definition of dt in fieldfft_meep.m.
+
+;; Run it:
+(define print-field (lambda () (print "field:, " (meep-time) ", "
+	       (real-part (get-field-point Ex (vector3 0 0 (+ (* -0.5 sz) Lcav (* 0.5 dpad))))) "\n")))
+(run-until endt (after-time (- endt 250) print-field))

scheme/meep.scm.in

@@ -124,19 +124,20 @@
 ; ****************************************************************
 ; Multilevel-atom nonlinear susceptibilities
 
-;(define-class transition no-parent
-;  (define-property from-level no-default 'integer non-negative?)
-;  (define-property to-level no-default 'integer non-negative?)
-;  (define-property transition-rate 0 'number) ; nonradiative rate (0 if none)
-;  (define-property frequency 0 'number) ; radiative frequency (0 if none)
-;  (define-property sigma-diag (vector3 1 1 1) 'vector3) ; per-transition sigma
-;  (define-property gamma 0 'number)) ; optical damping rate
-
-;(define (transition-time t) (transition-rate (/ t)))
-
-;(define-class multilevel-atom susceptibility
-;  (define-property initial-populations '() (make-list-type 'number))
-;  (define-property transitions '() (make-list-type 'transition)))
+(define-class transition no-parent
+  (define-property from-level no-default 'integer non-negative?)
+  (define-property to-level no-default 'integer non-negative?)
+  (define-property transition-rate 0 'number) ; nonradiative rate (0 if none)
+  (define-property frequency 0 'number) ; radiative frequency (0 if none)
+  (define-property sigma-diag (vector3 1 1 1) 'vector3) ; per-transition sigma
+  (define-property gamma 0 'number) ; optical damping rate
+  (define-property pumping-rate 0 'number)) ; pumping rate (0 if none)
+
+(define (transition-time t) (transition-rate (/ t)))
+
+(define-class multilevel-atom susceptibility
+  (define-property initial-populations '() (make-list-type 'number))
+  (define-property transitions '() (make-list-type 'transition)))
 
 ; ****************************************************************
 ; Add some predefined variables, for convenience:

scheme/structure.cpp

@@ -1148,11 +1148,9 @@ static bool susceptibility_equiv(const susceptibility *o0,
 				 const susceptibility *o)
 {
 if (o0->which_subclass != o->which_subclass) return 0;
-#if 0
 if (o0->which_subclass == susceptibility::MULTILEVEL_ATOM) {
 if (!multilevel_atom_equal(o0->subclass.multilevel_atom_data, o->subclass.multilevel_atom_data)) return 0;
 }
-#endif
 else if (o0->which_subclass == susceptibility::DRUDE_SUSCEPTIBILITY) {
 if (!drude_susceptibility_equal(o0->subclass.drude_susceptibility_data, o->subclass.drude_susceptibility_data)) return 0;
 }
@@ -1215,7 +1213,6 @@ void geom_epsilon::sigma_row(meep::component c, double sigrow[3],
     material_type_destroy(material);
 }
 
-#if 0
 /* make multilevel_susceptibility from scheme input data */
 static meep::susceptibility *make_multilevel_sus(const multilevel_atom *d) {
   if (!d || d->transitions.num_items == 0) return NULL;
@@ -1249,8 +1246,8 @@ static meep::susceptibility *make_multilevel_sus(const multilevel_atom *d) {
   for (int t = 0; t < d->transitions.num_items; ++t) {
     int i = d->transitions.items[t].from_level - minlev;
     int j = d->transitions.items[t].to_level - minlev;
-    Gamma[i*L+i] += d->transitions.items[t].transition_rate;
-    Gamma[j*L+i] -= d->transitions.items[t].transition_rate;
+    Gamma[i*L+i] += + d->transitions.items[t].transition_rate + d->transitions.items[t].pumping_rate;
+    Gamma[j*L+i] -= + d->transitions.items[t].transition_rate + d->transitions.items[t].pumping_rate;
   }
 
   // initial populations of each level
@@ -1264,6 +1261,9 @@ static meep::susceptibility *make_multilevel_sus(const multilevel_atom *d) {
   meep::realnum *omega = new meep::realnum[T];
   meep::realnum *gamma = new meep::realnum[T];
   meep::realnum *sigmat = new meep::realnum[T * 5];
+
+  const double pi = 3.14159265358979323846264338327950288; // need pi below.
+
   for (int t = 0, tr = 0; t < d->transitions.num_items; ++t)
     if (d->transitions.items[t].frequency != 0) {
       omega[tr] = d->transitions.items[t].frequency; // no 2*pi here
@@ -1280,8 +1280,8 @@ static meep::susceptibility *make_multilevel_sus(const multilevel_atom *d) {
       }
       int i = d->transitions.items[t].from_level - minlev;
       int j = d->transitions.items[t].to_level - minlev;
-      alpha[i * T + tr] = -1.0 / omega[tr];
-      alpha[j * T + tr] = +1.0 / omega[tr];
+      alpha[i * T + tr] = -1.0 / (2*pi*omega[tr]); // but we *do* need the 2*pi here. -- AWC
+      alpha[j * T + tr] = +1.0 / (2*pi*omega[tr]); 
       ++tr;
     }
 
@@ -1297,7 +1297,6 @@ static meep::susceptibility *make_multilevel_sus(const multilevel_atom *d) {
 
   return s;
 }
-#endif
 
 // add a polarization to the list if it is not already there
 static pol *add_pol(pol *pols, const susceptibility *user_s)
@@ -1384,13 +1383,11 @@ void geom_epsilon::add_susceptibilities(meep::field_type ft,
 	}
 	break;
       }
-#if 0
       case susceptibility::MULTILEVEL_ATOM: {
 	multilevel_atom *d = p->user_s.subclass.multilevel_atom_data;
 	sus = make_multilevel_sus(d);
 	break;
       }
-#endif
       default:
 	meep::abort("unknown susceptibility type");
     }

src/multilevel-atom.cpp

@@ -23,8 +23,6 @@
 #include "meep_internals.hpp"
 #include "config.h"
 
-using namespace std;
-
 namespace meep {
 
 multilevel_susceptibility::multilevel_susceptibility(int theL, int theT,
@@ -259,9 +257,10 @@ void multilevel_susceptibility::update_P
 
     // Ntmp = (I - Gamma * dt/2) * N
     for (int l1 = 0; l1 < L; ++l1) {
-      Ntmp[l1] = (1.0 - Gamma[l1*L + l1]*dt2) * N[l1]; // diagonal term
-      for (int l2 = 0; l2 < l1; ++l2) Ntmp[l1] -= Gamma[l1*L+l2]*dt2 * N[l2];
-      for (int l2 = l1+1; l2 < L; ++l2) Ntmp[l1] -= Gamma[l1*L+l2]*dt2 * N[l2];
+      Ntmp[l1] = 0;
+      for (int l2 = 0; l2 < L; ++l2) {
+	Ntmp[l1] += ((l1 == l2) - Gamma[l1*L+l2]*dt2) * N[l2];
+      }
     }
 
     // compute E*8 at point i
@@ -281,22 +280,31 @@ void multilevel_susceptibility::update_P
 
     // Ntmp = Ntmp + alpha * E * dP
     for (int t = 0; t < T; ++t) {
-      // compute 32 * E * dP at point i
+      // compute 32 * E * dP and 64 * E * P at point i
       double EdP32 = 0;
+      double EPave64 = 0;
+      double gperpdt = gamma[t]*pi*dt;
       for (idot = 0; idot < 3 && cdot[idot] != Dielectric; ++idot) {
 	realnum *p = d->P[cdot[idot]][0][t], *pp = d->P_prev[cdot[idot]][0][t];
 	realnum dP = p[i]+p[i+o1[idot]]+p[i+o2[idot]]+p[i+o1[idot]+o2[idot]]
 	  - (pp[i]+pp[i+o1[idot]]+pp[i+o2[idot]]+pp[i+o1[idot]+o2[idot]]);
+	realnum Pave2 = p[i]+p[i+o1[idot]]+p[i+o2[idot]]+p[i+o1[idot]+o2[idot]]
+	  + (pp[i]+pp[i+o1[idot]]+pp[i+o2[idot]]+pp[i+o1[idot]+o2[idot]]);
 	EdP32 += dP * E8[idot][0];
+	EPave64 += Pave2 * E8[idot][0];
 	if (d->P[cdot[idot]][1]) {
 	  p = d->P[cdot[idot]][1][t]; pp = d->P_prev[cdot[idot]][1][t];
 	  dP = p[i]+p[i+o1[idot]]+p[i+o2[idot]]+p[i+o1[idot]+o2[idot]]
+	    - (pp[i]+pp[i+o1[idot]]+pp[i+o2[idot]]+pp[i+o1[idot]+o2[idot]]);
+	  Pave2 = p[i]+p[i+o1[idot]]+p[i+o2[idot]]+p[i+o1[idot]+o2[idot]]
 	    + (pp[i]+pp[i+o1[idot]]+pp[i+o2[idot]]+pp[i+o1[idot]+o2[idot]]);
 	  EdP32 += dP * E8[idot][1];
+	  EPave64 += Pave2 * E8[idot][1];
 	}
       }
       EdP32 *= 0.03125; /* divide by 32 */
-      for (int l = 0; l < L; ++l) Ntmp[l] += alpha[l*T + t] * EdP32;
+      EPave64 *= 0.015625; /* divide by 64 (extra factor of 1/2 is from P_current + P_previous) */ 
+      for (int l = 0; l < L; ++l) Ntmp[l] += alpha[l*T+t]*EdP32 + alpha[l*T+t]*gperpdt*EPave64;
     }
 
     // N = GammaInv * Ntmp
@@ -308,10 +316,12 @@ void multilevel_susceptibility::update_P
 
   // each P is updated as a damped harmonic oscillator
   for (int t = 0; t < T; ++t) {
-    const double omega2pi = 2*pi*omega[t], g2pi = gamma[t]*2*pi;
-    const double omega0dtsqr = omega2pi * omega2pi * dt * dt;
-    const double gamma1inv = 1 / (1 + g2pi*dt/2), gamma1 = (1 - g2pi*dt/2);
-
+    const double omega2pi = 2*pi*omega[t], g2pi = gamma[t]*2*pi, gperp = gamma[t]*pi;
+    const double omega0dtsqrCorrected = omega2pi*omega2pi*dt*dt + gperp*gperp*dt*dt;
+    const double gamma1inv = 1 / (1 + g2pi*dt2), gamma1 = (1 - g2pi*dt2);
+    const double dtsqr = dt*dt;
+    // note that gamma[t]*2*pi = 2*gamma_perp as one would usually write it in SALT. -- AWC
+
     // figure out which levels this transition couples
     int lp = -1, lm = -1;
     for (int l = 0; l < L; ++l) {
@@ -350,11 +360,11 @@ void multilevel_susceptibility::update_P
 	    realnum pcur = p[i];
 	    const realnum *Ni = N + i*L;
 	    // dNi is population inversion for this transition
-	    double dNi = -0.25 * (Ni[lp]+Ni[lp+o1]+Ni[lp+o2]+Ni[lp+o1+o2]
-				  -Ni[lm]-Ni[lm+o1]-Ni[lm+o2]-Ni[lm+o1+o2]);
-	    p[i] = gamma1inv * (pcur * (2 - omega0dtsqr)
+	    double dNi = 0.25 * (Ni[lp]+Ni[lp+o1]+Ni[lp+o2]+Ni[lp+o1+o2]
+				 -Ni[lm]-Ni[lm+o1]-Ni[lm+o2]-Ni[lm+o1+o2]);
+	    p[i] = gamma1inv * (pcur * (2 - omega0dtsqrCorrected) 
 				- gamma1 * pp[i]
-				+ omega0dtsqr * (st * s[i] * w[i])) * dNi;
+				- dtsqr * (st * s[i] * w[i]) * dNi);
 	    pp[i] = pcur;
 	  }
 	}