Field normalization of first band at zero wave vector

[thchr]

Hi Steven,

When checking the normalization of E and H I get the following odd behavior when k-points is (vector3 0 0 0) . Using the .ctl-file below, I calculate the overlap between the nth and mth bands (in the D_n^†E_m and B_n^†H_m sense) and get (TM-calculation with µ = 1 and ε = 13 cylinders; only real part of overlap shown):

---DE overlaps---
             m = 1,       m = 2,       m = 3,       m = 4,       m = 5
n = 1 |   0.00    ,    0.00    ,    0.00    ,    0.00    ,    0.00    
n = 2 |   0.00    ,     1.0    ,    2.64E-16,    2.01E-16,    2.91E-16
n = 3 |   0.00    ,    2.36E-16,     1.0    ,    1.39E-17,   -2.78E-17
n = 4 |   0.00    ,    2.01E-16,    1.39E-17,     1.0    ,    3.38E-16
n = 5 |   0.00    ,    3.19E-16,   -3.47E-17,    3.38E-16,     1.0    

---BH overlaps---
             m = 1,       m = 2,       m = 3,       m = 4,       m = 5
n = 1 |    1.0    ,    2.51E-19,    9.02E-19,    2.68E-19,   -2.51E-18
n = 2 |   2.51E-19,    1.00    ,    5.00E-16,   -6.83E-16,    8.90E-16
n = 3 |   9.02E-19,    5.00E-16,    1.00    ,    1.39E-15,   -4.20E-16
n = 4 |   2.68E-19,   -6.83E-16,    1.39E-15,     1.0    ,    5.00E-16
n = 5 |  -2.51E-18,    8.90E-16,   -4.20E-16,    5.00E-16,    1.00    

The issue is in the left column/top row of the D^†E calculations. For more complicated geometries, I also get strange results for the B^†H term.
[If the permeability mu is not a unity matrix, this column/row is not zero, but a small value (10^-2-10^-6); in complicated geometries the D^†E+B^†H overlap sums to ~6 instead of 2]

If I change k-points to, say, (vector3 1e-4 0 0) the problem is remedied:

---DE overlaps---
             m = 1,       m = 2,       m = 3,       m = 4,       m = 5
n = 1 |   1.00    ,   -1.92E-12,    3.91E-12,    2.00E-12,    9.02E-13
n = 2 |  -1.92E-12,    1.00    ,    1.11E-16,   -3.61E-16,    1.73E-16
n = 3 |   3.91E-12,    1.11E-16,    1.00    ,   -4.72E-16,   -1.67E-16
n = 4 |   2.00E-12,   -3.61E-16,   -4.72E-16,    1.00    ,    2.01E-16
n = 5 |   9.02E-13,    1.73E-16,   -1.11E-16,    2.01E-16,     1.0    

---BH overlaps---
             m = 1,       m = 2,       m = 3,       m = 4,       m = 5
n = 1 |    1.0    ,    4.89E-12,    3.69E-11,   -2.37E-11,    1.84E-10
n = 2 |   4.89E-12,     1.0    ,   -1.39E-09,    3.99E-09,   -6.80E-10    
n = 3 |   3.69E-11,   -1.39E-09,     1.0    ,   -2.23E-09,    2.12E-11   
n = 4 |  -2.37E-11,    3.99E-09,   -2.23E-09,     1.0    ,   -2.63E-09   
n = 5 |   1.84E-10,   -6.80E-10,    2.12E-11,   -2.63E-09,     1.0      

Is this issue arising because the solution to the k = 0 case is trivial (ω = 0)? I suppose it is not a big issue, but it is slightly worrying e.g. for super cell calculations if overlap integrals between all bands are needed (= what I am doing).

Thanks!
-Thomas

.ctl file

(the B^†H overlap is left out; trivial repetitions)

(use-modules (ice-9 format)) ; load 'format' to allow precise control of print output

; define default variable values editable from command-line call
(define-param timesym "timebroken")

(set! num-bands 5)
(set! resolution 32)

; define material properties [based on input timesym, evaluated by if-statment (cond)]
(define matprops
(cond ((string=? timesym "timeinvariant")  ; time-invariant isotropic dielectric response
          (make dielectric (epsilon 13)))
      ((string=? timesym "timebroken")     ; time-broken anisotropic response
          (make medium-anisotropic (epsilon-diag 13 13 13) (epsilon-offdiag 0 0 0)
                                                    (mu-diag 1 1 1) (mu-offdiag 0+0.40i 0 0) )
      )
))

; define a simple cylinder-in-square-lattice geometry
(set! geometry-lattice (make lattice (size 1 1 no-size)   (basis1 1 0) (basis2 0 1)))
(set! geometry (list (make cylinder (center 0 0 0) (radius 0.1) (height infinity)
                                     (material matprops))))
 
(set! k-points (list (vector3 0 0 0))) ; k-point at Gamma

; function to compute band overlaps
(define (get-overlaps)
   ; function to calculate overlap between D_n and E_m or B_n and H_m
   (define (DEorBH-integrand r DB-n EH-m)  (vector3-cdot DB-n EH-m) )
   ; create a list for loop (1,2,3,...,num-bands)
   (define iter-list (list-tail (iota (+ 1 num-bands)) 1))

   ; loop over the list twice and calculate the required integrals along the way; first DE terms
   (print "\n\n---DE overlaps---\n\t  ")
   (map (lambda (m) (print "   m = " m ",    " )) iter-list)
   (map  (lambda (n)
            (print "\nn = " n " |\t")
            (get-dfield n)
            (let  ((D-n (field-copy cur-field)))
                  (map  (lambda (m)
                           (get-efield m)
                           (let  ((E-m (field-copy cur-field)))
                                 (let  ((overlap (integrate-fields DEorBH-integrand D-n E-m)))
                                       (print (format #f "~10,2,2g" (real-part overlap)) ",  ")
                                 )
                           )
                        )
                  iter-list )
            )
         )
    iter-list )
   (print "\n\n") ; clear a line
)

(run-tm get-overlaps) ; run tm-calculation, calculate overlaps at the end

Separate normalization question

Perhaps I should be opening a separate issue for this, but it's more of a question than an issue: I was wondering if the normalization-choices: ∫ D^†E dr = ∫ B^†H dr = 1 imply that E and H are normalized independently of each other? I.e. are ∇×E ≠ iωB and ∇×H ≠ -iωD?
I'm interested in calculating integrals of the type

<n|f|m> = ∫ F_n^†(r)M(r)f(r)F_m(r) dr

with the 6-vector F = (E; H), some scalar function f, and the usual weight matrix M = diag(ε, μ). If E and H are normalized independently, it seems to me that a naive construction of F via get-efield and get-hfield would lead to trouble. Is this correctly understood?

I don't suppose there's an option which would enable a normalization in the 6-vector sense, i.e. in a sense such that ε|E|² + µ|H|² integrates to unity, which is simultaneously consistent with ∇×E = iωB and ∇×H = -iωD?

[stevengj]

The basic issue here is that at k=0 the solution for the first band (the "electrostatic" solution at ω=0) is ambiguous. The E and H fields are completely decoupled there, and the choice of solution is therefore somewhat arbitrary.

You can get a unique quasistatic solution by taking the limit as k goes to 0 from some direction. (MPB can't do this for you, because it depends on which direction you approach k=0 from.)