Differential Boundary Conditions in Electrodynamics

Terje Vold

Continuum Technology, Inc.

The Problem

Use of Geometric Algebra and Calculus

The Usual Solution

A Solution using Differential Boundary Conditions

Plans

The Problem

Numerically simulate electromagnetic fields 
in general situations of interest that are well suited to surface integral methods 
with the best accuracy and least computational effort that is practical.

Characteristics of these types of problems 
	space is filled with matter having linear electric and magnetic polarizabilities
	the matter is piecewise homogeneous (lumps or pieces of homogeneous material)
	the surface to volume ratio is not large
	the matter is stationary though generally has changing electric charge and current

Examples
	static and low frequency magnetic shielding
	high frequency electromagnetic shielding
	optical integrated circuits
	antennas

Why is this a problem?   
	Implementation details are complicated, and
	Faster and more accurate variations of the general solution method are always good

The solution methods considered here are usually called surface integral methods and use Green's functions.  The general idea is:

a) If we have the correct field values on the boundary of a physically homogeneous region of space, 
then we can calculate the field at any point inside as an integral over the boundary or "surface".

b) The field values on the 2D boundaries are determined by requiring them to satisfy boundary conditions (and also meet self-consistency requirements).

In this way, the original problem of determining field values throughout 3D space is reduced to determining the correct field values only on 2D boundaries.

Use of Geometric Algebra and Calculus

This work uses geometric algebra, a very useful algebra that includes conventional vector algebra as a subsystem.

Quick Introduction to Geometric Algebra

Analogy using complex algebra

In rough terms, one might say that 	conventional vector algebra is to geometric algebra
				as 	real algebra is to complex algebra.  
				
For example,  

Real:  		Cos[θ+ϕ] = Cos[θ]Cos[ϕ] - Sin[θ]Sin[ϕ]
		Sin[θ+ϕ] = Cos[θ]Sin[ϕ] + Sin[θ]Cos[ϕ]
		
Complex:  	above, plus     Exp[I (θ+ϕ)] = Exp[Iθ] Exp[iϕ]

Analogy using component, vector, and geometric algebra

Or one might say
	
	geometric algebra is to vector algebra 
	as vector algebra is to vector component algebra.
	
For example, Maxwell's equations in SI units:

a) as Maxwell would write using vector component algebra,

	(∂_x)(Subscript[D, x])+(∂_x)(Subscript[D, x])+(∂_x)(Subscript[D, x]) = ρ

	(∂_z)(Subscript[H, y])-(∂_y)(Subscript[H, z])-(∂_t)(Subscript[D, x]) = (Subscript[J, x])
	(∂_x)(Subscript[H, z])-(∂_z)(Subscript[H, x])-(∂_t)(Subscript[D, y]) = (Subscript[J, y])
	(∂_y)(Subscript[H, x])-(∂_x)(Subscript[H, y])-(∂_t)(Subscript[D, z]) = (Subscript[J, z])

	(∂_z)(Subscript[E, y])-(∂_y)(Subscript[E, z])+(∂_t)(Subscript[B, x]) = 0
	(∂_x)(Subscript[E, z])-(∂_z)(Subscript[E, x])+(∂_t)(Subscript[B, y]) = 0
	(∂_y)(Subscript[E, x])-(∂_x)(Subscript[E, y])+(∂_t)(Subscript[B, z]) = 0

	(∂_x)(Subscript[B, x])+(∂_x)(Subscript[B, x])+(∂_x)(Subscript[B, x]) = 0

b) as W.J.Gibbs would write using his vector algebra, 
		(OverVector[∇])·(OverVector[D]) = ρ
		(OverVector[∇])×(OverVector[H]) - (∂_t)(OverVector[D]) = (OverVector[J])
		(OverVector[∇])×(OverVector[E]) + (∂_t)(OverVector[B]) = 0
		(OverVector[∇])·(OverVector[B]) = 0

c) as Hestenes writes using geometric algebra, 
		(OverVector[∇])G + (1/c)(∂_t)G = ρ-(OverVector[J])

where  G = (OverVector[D]) + III (OverVector[H])/c.

Geometric Algebra of 3D space in a nutshell

The basic ("atomic") elements of the algebra are 
scalars (real or complex numbers), vectors, bivectors (vectors times III), trivectors (scalars times III).   
III is the "unit trivector".   It satisfies  (III^2)=-1, but is distinct from the imaginary unit i.  

Any elements of the algebra can be added or multiplied together,   
just as you can add or multiply real and imaginary numbers. 
Such a general element is called a "multivector."

Addition is unsurprising. 

Multiplication is associative but NOT commutative.  
The even or symmetric part of multiplication between two vectors is the dot product. 
The odd or antisymmetric part of multiplication between two vectors is the "wedge product";
	this wedge product times III equals the cross product.

Examples of use of Geometric Algebra and Calculus

Introduction

Lots of subtleties and difficulties with symbolic calculation in geometric algebra --- 
a) implementing a noncommutative product 
b) finding notation that is acceptably consistent with tradition and works simply with Mathematica
c) implementing the noncommutative vector derivative - which might operate to the left or right
d) deciding what transformation rules to always apply and what to apply only when specified

A frame-dependent numeric implementation 
using 8x8x8 arrays and 8-element lists of multivector components 
is straightforward and convenient.

NotebookDirectory[]

NotebookDirectory[]<>"Geometric Algebra.m"

Get[NotebookDirectory[]<>"Geometric Algebra.m"]

Dipole Field example

Static dipole:

staticDipolePotential=(((OverVector[m])·(OverVector[r]) )/((((OverVector[r])·(OverVector[r])))^(3/2)))

(OverVector[▽])∘staticDipolePotential

Exact time-varying dipole:

greensfunction[k_,R_]:=(1/(4π ((\[LeftBracketingBar]R\[RightBracketingBar])^3)))Exp[I k \[LeftBracketingBar]R\[RightBracketingBar]](R+I k \[LeftBracketingBar]R\[RightBracketingBar](\[LeftBracketingBar]R\[RightBracketingBar]-R))

DipoleField[k_,d_?Grade1Q,r_?Grade1Q]:=(d·(Subscript[(OverVector[D]), r]))⋄greensfunction[k,-r]-I k  greensfunction[k,-r]∘d

DipoleField[k,(OverVector[d]),(OverVector[r])]

(DipoleField[k,(OverVector[d]),(OverVector[r])]//ExpandGP)/.ToCrossRule

Grade[k]={0};

(DipoleField[k,(OverVector[d]),(OverVector[r])]//ExpandGP)/.ToCrossRule

Directional Derivative example

((OverVector[n])·(Subscript[(OverVector[▽]), (OverVector[r])]))⋄DipoleField[k,(OverVector[d]),(OverVector[r])]

%/.ToCrossRule

Matrix giving directional derivative as a linear function of dipole moment components {dx,dy,dz} ---

DDDF=(((OverVector[n])·(Subscript[(OverVector[▽]), (OverVector[r])]))⋄DipoleField[k,(OverVector[d]),(OverVector[r])]/.ToCrossRule)/.{(OverVector[d])->{dx,dy,dz},(OverVector[n])->{nx,ny,nz},(OverVector[r])->{rx,ry,rz}}

CoefficientList[DDDF,III]//Dimensions

fieldComponentsList=Flatten[Transpose[CoefficientList[DDDF,III]]];
fieldComponentsList//Dimensions
fieldComponentsList//MatrixForm

FieldFromMomentMatrix=CoefficientArrays[fieldComponentsList,{dx,dy,dz}][[2]];
FieldFromMomentMatrix//Dimensions
FieldFromMomentMatrix//MatrixForm

Maxwell's Equations example

Grade[k]={0};
Grade[ρ]={0};

Clear[F]

Clear[OverVector]

FRule={F->(OverVector[ℰ])[(OverVector[r])]+III (OverVector[ℬ])[(OverVector[r])]};

maxwell=((OverVector[▽])∘F-I k F  == ρ-(OverVector[J]))

maxwell/.FRule

maxwellExpanded=ExpandGP[maxwell/.FRule]

maxwellExpandedUsingCross=maxwellExpanded/.ToCrossRule

GradeProjection[maxwellExpandedUsingCross,0]

GradeProjection[maxwellExpandedUsingCross,1]

GradeProjection[maxwellExpandedUsingCross,2]

GradeProjection[maxwellExpandedUsingCross,3]

OverVector/:(OverVector[▽])∘a__∘n_.:=VD[a∘n,DefaultDerivativeArgument]

The Usual Solution

Green's Integral

The electromagnetic field G = (OverVector[D]) + III(OverVector[H])/c  inside any bounded volume is given by the field G on the surface and any enclosed charge and current sources S = ρ - (OverVector[J]) by Green's integral,

G[(OverVector[r])']=(∫_(∂V))g[(OverVector[r])-(OverVector[r])'] d(OverVector[a]) G[(OverVector[r])] - (∫_V)g[(OverVector[r])-(OverVector[r])'] da S[(OverVector[r])]

where

g[(OverVector[R])]=(1/(4π((|(OverVector[R])|)^3)))Exp[i k |(OverVector[R])|]((OverVector[R])+i k|(OverVector[R])|(|(OverVector[R])|-(OverVector[R]) ))

It is often convenient to consider the second integral to be an "incident" or "impressed" field (G^i).  With this, and explicitly writing G as the electric vector field plus the magnetic bivector field, and then defining "effective" surface charges and currents, we have

G[(OverVector[r])'] = (∫_(∂V))g[(OverVector[r])-(OverVector[r])'] d(OverVector[a]) G[(OverVector[r])]  +  (G^i)[(OverVector[r])']

	=(∫_(∂V))g[(OverVector[r])-(OverVector[r])'](d(OverVector[a]) ·((OverVector[D])[(OverVector[r])]+III (OverVector[H])[(OverVector[r])]) + d(OverVector[a]) ⋀((OverVector[D])[(OverVector[r])]+III (OverVector[H])[(OverVector[r])]) ) + (G^i)[(OverVector[r])']

	=(∫_(∂V))g[(OverVector[r])-(OverVector[r])'](d(OverVector[a]) ·(OverVector[D])[(OverVector[r])]-d(OverVector[a])× (OverVector[H])[(OverVector[r])] +III d(OverVector[a])×(OverVector[D])[(OverVector[r])]+III d(OverVector[a])· (OverVector[H])[(OverVector[r])]) + (G^i)[(OverVector[r])']

	= (∫_(∂V))g[(OverVector[r])-(OverVector[r])']da ((Subscript[σ, e])[(OverVector[r])]- (Subscript[(OverVector[κ]), e])[(OverVector[r])] -III (Subscript[(OverVector[κ]), m])[(OverVector[r])]+III (Subscript[σ, m])[(OverVector[r])]) + (G^i)[(OverVector[r])']

where the effective surface charges and currents (Subscript[σ, e]), (Subscript[(OverVector[κ]), e]), (Subscript[σ, m]), (Subscript[(OverVector[κ]), m])  are just suggestive shorthand for the corresponding quantities in the preceeding line.

This says that G at any point r' can be calculated as the applied field plus an integral over the boundary, where the integrand has terms representing effective electric charge, effective electric current, effective magnetic current, and effective magnetic charge.

Boundary Conditions

The usual method is to find effective surface charges and currents that give rise to fields that satisfy the well-known boundary conditions between two regions, say A and B, that follow from Maxwell's equations:

(OverHat[n])·(Subscript[(OverVector[D]), A])- (OverHat[n])·(Subscript[(OverVector[D]), B]) = 0
 (OverHat[n])×(Subscript[(OverVector[H]), A])-  (OverHat[n])×(Subscript[(OverVector[H]), B]) = 0
	 (OverHat[n])×(Subscript[(OverVector[E]), A])- (OverHat[n])×(Subscript[(OverVector[E]), B]) = 0
 (OverHat[n])·(Subscript[(OverVector[B]), A])-  (OverHat[n])·(Subscript[(OverVector[B]), B]) = 0

The simplest method of solution at this point is to parameterize the boundary fields (or equivalent effective boundary charge and current) using some reasonable set of "basis" functions, and then choose the parameter values by requiring that the Maxwell boundary conditions are satified at some set of discrete points on the boundary.

A more complex variation is to require the boundary conditions by integrating the point conditions over the boundary, weighted by the elements of a set of typically localized, peaked weighting functions.  This may be significantly more complex, and does not always give better results.  This (or a particular variant) is called the "method of moments" in electrical engineering.

In either case, the result is a matrix equation to be solved for some finite set of parameters.

A Solution using Differential Boundary Conditions

Solution Method Sketch

Rather that using boundary conditions at discrete points, or using "average" boundary conditions obtained by integrating the boundary conditions over small subsections of the boundary, we consider requiring 

	a) the usual point boundary condition
	    plus a tangential derivative boundary condition, 
	b) field parameterization by point charge and current sources (dipoles) 
	    that are _not_ on the boundary.

This has advantages: 
	a) is simpler than integrals over the boundary condition 
	b) is simpler than integrating over surface charge and current

But this also has disadvantages:
	a) more computation for the tangential derivative boundary condition
	b) awkward handling of thin homogeneous volumes due to the field parameterization method

Boundary Condition for Field Derivatives for Curved Boundaries

The usual boundary conditions between dissimilar polarizable media A and C with boundary normal (OverHat[n]) may be written

(OverHat[n]) (OverHat[n])⋀((1/(Subscript[η, A]))(Subscript[G, A])-(1/(Subscript[η, C]))(Subscript[G, C]))  +  (OverHat[n])  (OverHat[n])·((Subscript[G, A])-(Subscript[G, C]))=0

where

(Subscript[G, A])=(Subscript[(OverVector[D]), A])+III (Subscript[(OverVector[H]), A])

(Subscript[F, A])=(Subscript[(OverVector[E]), A])+III (Subscript[(OverVector[B]), A])

(Subscript[G, A])=(Subscript[η, A])(Subscript[F, A])

A general boundary condition for field derivatives, with zero boundary charge and current, is

(OverHat[n])⋀(OverVector[∇])((OverHat[n]) (OverHat[n])⋀((1/(Subscript[η, A]))(Subscript[G, A])-(1/(Subscript[η, C]))(Subscript[G, C]))  +  (OverHat[n])  (OverHat[n])·((Subscript[G, A])-(Subscript[G, C])))=0

For a curved boundary, n depends on position.  We express this as follows.

Let { (Subscript[(OverHat[e]), 1]), (Subscript[(OverHat[e]), 2]), (Subscript[(OverHat[e]), 3])=(OverHat[n ]) } be an orthonormal frame, with {(Subscript[(OverHat[e]), 1]), (Subscript[(OverHat[e]), 2])} chosen to be the curvature unit eigenvectors associated with (diagonal) curvature matrix c.  Then

(OverVector[r]) ≈(OverVector[ρ]) = (Superscript[ρ, i])(Subscript[(OverHat[e]), i])

(OverHat[n]) ≈ (Subscript[(OverHat[e]), 3])+ (Superscript[ρ, i])(Subscript[c, i]^j)(Subscript[(OverHat[e]), j])

(Subscript[G, X])[(OverVector[r])]≈ (Subscript[G, X])+(OverVector[ρ])·(OverVector[∇])(Subscript[G, X]) = (Subscript[G, X])+(Superscript[ρ, i])(∂_i)(Subscript[G, X])

and

(OverHat[n])⋀(OverVector[∇]) = (Subscript[(OverHat[e]), 3])((OverHat[e])^i)(∂_i)

where sums over {i,1,2} and not i=3 are implied, and cij is a 2x2 curvature matrix that is diagonal in the eigenvector frame.

We can show that the general boundary condition for field derivatives written above may be expanded as

-((OverHat[n])  (OverHat[n])·((1/(Subscript[η, A]))(OverHat[n])·(OverVector[∇]) (Subscript[G, A]))  +  (OverHat[n])  (OverHat[n])⋀((OverHat[n])·(OverVector[∇]) (Subscript[G, A]))) - (1-(1/(Subscript[η, A]))) 𝒞[(Subscript[G, A])]  =  -((OverHat[n])  (OverHat[n])·((1/(Subscript[η, C]))(OverHat[n])·(OverVector[∇]) (Subscript[G, C]))  +  (OverHat[n])  (OverHat[n])⋀((OverHat[n])·(OverVector[∇]) (Subscript[G, C]))) - (1-(1/(Subscript[η, C]))) 𝒞[(Subscript[G, C])]

where C is the traceless symmetric linear operator

𝒞[M]= (Superscript[(OverHat[e]), i])(Subscript[c, i]^j) (Subscript[(OverHat[e]), j])·M - Tr[c](Subscript[(OverHat[e]), 3]) (Subscript[(OverHat[e]), 3])·M

A goal is to have Mathematica do these symbolic manipulations for me!   I'm working on it.

Some preliminary numerical results

Sphere of dielectric represented by 6 boundary points

A sphere of dielectric material, represented by six boundary points forming an octahedron, 
with an applied (or "incident" or "impressed") electromagnetic wave.

grData1=First[PolyhedronData["Octahedron"]]

Graphics3D[{Opacity[0.3],grData1,Opacity[0.7],Sphere[{0,0,0},1/Sqrt[2]]}]

dipolePoints={{{-0.10455179533282444`,0.`,0.`},{-4.4028315226269905`,0.`,0.`}},{{0.`,0.10455179533282444`,0.`},{0.`,4.4028315226269905`,0.`}},{{0.`,0.`,-0.10455179533282444`},{0.`,0.`,-4.4028315226269905`}},{{0.`,0.`,0.10455179533282444`},{0.`,0.`,4.4028315226269905`}},{{0.`,-0.10455179533282444`,0.`},{0.`,-4.4028315226269905`,0.`}},{{0.10455179533282444`,0.`,0.`},{4.4028315226269905`,0.`,0.`}}};

Graphics3D[{Point[Flatten[dipolePoints,1]],Opacity[0.3],grData1,Opacity[0.7],Sphere[{0,0,0},1/Sqrt[2]]}]

Field with intermediate frequency incident wave, 
permittivity ϵ = 3

Show[plot2,plot34,Graphics3D[{Opacity[0.2],Sphere[{0,0,0},TetrahedronRadius]}]]

Dielectric Sphere represented by 6 points, freq=0, permittivity ϵ = 100

Exact theoretical result:  Ein = (3/(2+ϵ))Eapp,  Eoutpole =  ϵ Ein. 
So for ϵ=100 and Eapp=1,  expect 
	Ein = 3/102 and so Din = 300/102≈2.94;  
	Eoutpole = Doutpole = Din; 
	Eoutequator=0.0294

Partition[Partition[Partition[GofMMatrix[].mainsolu,6],2],2]//Chop//MatrixForm

Chop[(({{(({{0, 0, 1-0.9708840501119995`, 0, 0, 0}, {0, 0, 4.151574978276672`, 0, 0, 0}})), (({{0, 0, 2.9115949888035777`, 0, 0, 0}, {0, 0, 0.07513402818526049`, 0, 0, 0}}))}, {(({{0, 0, 1-0.9708840501120541`, 0, 0, 0}, {0, 0, 4.151574978277289`, 0, 0, 0}})), (({{0, 0, 2.9115949888037664`, 0, 0, 0}, {0, 0, 0.07513402818594453`, 0, 0, 0}}))}, {(({{0, 0, 1+1.9261949606156785`, 0, 0, 0}, {0, 0, -8.19234765979352`, 0, 0, 0}})), (({{0, 0, 2.92619496061537`, 0, 0, 0}, {0, 0, 0.1416247118511897`, 0, 0, 0}}))}, {(({{0, 0, 1+1.9261949606158015`, 0, 0, 0}, {0, 0, -8.19234765979421`, 0, 0, 0}})), (({{0, 0, 2.9261949606152955`, 0, 0, 0}, {0, 0, 0.1416247118510136`, 0, 0, 0}}))}, {(({{0, 0, 1-0.9708840501122107`, 0, 0, 0}, {0, 0, 4.151574978278085`, 0, 0, 0}})), (({{0, 0, 2.91159498880382`, 0, 0, 0}, {0, 0, 0.07513402818591802`, 0, 0, 0}}))}, {(({{0, 0, 1-0.9708840501121052`, 0, 0, 0}, {0, 0, 4.15157497827755`, 0, 0, 0}})), (({{0, 0, 2.9115949888037918`, 0, 0, 0}, {0, 0, 0.07513402818595738`, 0, 0, 0}}))}}))[[All,All,All,{1,2,3}]]]//MatrixForm[#,TableDepth->2]&

Each row in the table above gives 
{{exterior D field, exterior D field deriv}, {interior D field, interior D field deriv}}.  

Results are very good:  We expect   

	interiorfield/exteriorfield = 1/ϵ = 1/100 at "equatorial" points, 
	interiorfield/exteriorfield = 1  at "polar" points, 
	
and 	interior D field = 300/102 ≈ 2.94.

In comparison, using 6 points without differential boundary conditions 
gives errors for the interior D field of order  ~10~30%.

Dielectric Sphere, 6 points, freq=10^7, permittivity ϵ = 4

Field induced in response to an incident plane wave:

EMFieldPlot3D[Flatten[soluarray02,1],{100000,100000},30]

Icosahedron (12 boundary points), intermediate frequency

Plot of induced electric and magnetic fields for a dielectric sphere with permittivity=3
represented by 12 boundary points describing an icosahedron
in response to a low frequency incident plane wave

Show[plot2,plot34,Graphics3D[{Opacity[0.2],Sphere[{0,0,0},icosahedronRadius]}]]

A more complex object:  a compound lens

There are 3 volumes:  the upper "lens," the lower "lens," and the space outside of both.

Graphics3D[{ico4bot,Opacity[0.5],ico2mid,ico4top}]

Possible locations of effective dipole sources for the upper boundary -

Graphics3D[Point[Flatten[SpaceComplexData[{3,2},"DipolePoints",Space->spaceExample04],1]]]

Possible locations of effective dipole sources for the lower boundary -

Graphics3D[Point[Flatten[SpaceComplexData[{1,2},"DipolePoints",Space->spaceExample04],1]]]

Possible locations of effective dipole sources for the middle boundary -

Graphics3D[Point[Flatten[SpaceComplexData[{2,2},"DipolePoints",Space->spaceExample04],1]]]

Plans

Extend the geometric algebra package - now only for the geometric algebra of a 3D vector space - to include and smoothly articulate with the geometric algebra of 4D spacetime.

Expand applicability of the geometric derivative.

Include reciprocal frame vectors in the geometric algebra package.

Include rotors in the geometric algebra package.

Develop integration for the geometric algebra package.

Compare practicality and accuracy/computation time using 
i) only field value boundary conditions, and 
ii) both field value and derivative boundary conditions.

Better understand the relationship between 
minimizing the electrodynamic action  and  the method of moments 
for solving this type of problem.