Problems performing 1D-FDTD for dispersive dielectric uisng lorentz model in Julia

So I am trying to model a dispersive dielectric using the Lorentz model in Julia, more specifically I am trying to obtain the frequency-dependent reflectivity of the material. I am ultimately trying to replicate the results in this paper. I have the set of update equations:

[Update equations](

along with the coefficients:


I am trying to implement this system but the numerical result I obtain is different than the analytical one. For reference:

[Results from my current implementation](

For some context, my source is a Gaussian pulse with central frequency 26 THz and bandwidth 6 THz. I also use the parameters fmax, nmax, d (the thickness of the slab) to determine the minimum wavelength I want to be able to resolve with my cell size through the relations lambda=c/(fmax*nmax), dzi=lambda/20,N = ceil(d/dzi), dz = d/N. Now I obtain the result above when I put in fmax=10,nmax=1, and d=304.8 which are not the expected values. When I put in the actual values fmax=32 (as the maximum freqeuncy in the source), nmax=16 (since the analtical result gives n=16+16i near the resonance), d=38.1 I get the following:

[Wrong? result for the right values](

Now the solution for the actual values results in a sin^2 pattern just as in the case with no dispersion. The distance between the peaks for a linear material can be found analytically to be c/(2*nmax*d) and here the peculiar thing is that the distance between the peaks is the same as the mean of this expression for all the different frequencies (as nmax is a function of frequency now since the material is dispersive). I remain dumbfounded by these results and can’t figure out why my implementation is not working.

I have my code written in Julia below (note: I am also taking D/epsilon_0 Q/sqrt(epsilon_0) and J/sqrt(epsilon_0)):


	ps = 1e-12 #s
	μm = 1e-6 #m
	THz = 1e12 #Hz
	F_per_m = ps^-4*μm^3 #(ps⁴*A²)/(kg*μm³)
	H_per_m = ps^2/μm #(kg*μm)/(ps²*A²)

	#Grid Resolution (Wavelength)
	c = 299792458*ps/μm #μm/ps
	λₘᵢₙ = c/(nₘₐₓ*fₘₐₓ) #μm
	Nₗ = 20 #number of points per wavelength
	Δₗ = λₘᵢₙ/Nₗ #μm

	#Grid Resolution (Structure)
	Nₛ = 4
	Δₛ = d/Nₛ #μm
	#Initial Grid Resolution (Overall)
	Δzᵢ = min(Δₗ,Δₛ) #μm
	#Snap Grid to Critical Dimension(s)
	dᵣ = d #μm
	N = ceil(dᵣ/Δzᵢ)
	Δz = dᵣ/N #μm
	#Determine Size of the Grid
	Nᵣ = Int64(N+23) #number of total cells
	z = (1:Nᵣ)*Δz
	#Compute Position of Materials
	n₁ = 13 #first cell of the material
	n₂ = Int64(n₁+round(d/Δz)-1) #last cell of the material
	#Add Materials to the Grid
	μᵣ = 1 #relative permeability of material
	μₛ = 1 #relative permeability of source
	εₛ = 1 #relative permittivity of source
	nₛ = sqrt(μₛ*εₛ) #refractive index of source
	UR=fill(μₛ,Nᵣ) #relative permeability across the grid
	#Calculate Time Step
	nb = 1 #refractive index of the boundary
	Δt = (nb*Δz)/(2c) #ps
	τ = 1/(2*6) #1/(2fₘₐₓ) #ps
	t₀ = 6τ #ps
	tₚᵣₒₚ = (nₘₐₓ*Nᵣ*Δz)/c #ps
	T = 12τ+tc*tₚᵣₒₚ #ps
	Steps = Int64(ceil(T/Δt)) #number of total time steps
	#Compute the Source Functions
	ω = 26*2π #THz
	t = (0:Steps-1)*Δt #ps
	δt = (nₛ*Δz)/(2c)+Δt/2 #ps
	Eₛ = exp.(-((t.-t₀)./τ).^2).*exp.(im*ω.*(t.-t₀)) #located at source
	Hₛ = -sqrt(εₛ/μₛ).*exp.(-((t.-t₀.+δt)./τ).^2).*exp.(im*ω.*(t.-t₀.+δt)) #located at source-1
	#Parameters for Equation of Motion
	μ₀ = 4π*1e-7*H_per_m #H/m
	ϵ₀ = 1/(c^2*μ₀) #F/m
	ε₀, εₓ = ones(Nᵣ), ones(Nᵣ)
	Γ, Ω, α, β = zeros(Nᵣ), zeros(Nᵣ), zeros(Nᵣ), zeros(Nᵣ)
	ε₀[n₁:n₂].= 9.66
	εₓ[n₁:n₂].= 6.52
	Γ[n₁:n₂].= 0.2 #THz
	Ω[n₁:n₂].= 24.9 #THz
	Ω₀ = sqrt(ε₀[n₁]/εₓ[n₁]*Ω[n₁]^2) #THz
	#Compute Update Coefficients
	#UR *= nb/(2*ri) .*(sin.(π*fₘₐₓ*ri*Δz/c))./(sin.(π*fₘₐₓ*Δt))
	#ε₀ *= nb/(2*ri) .*(sin.(π*fₘₐₓ*ri*Δz/c))./(sin.(π*fₘₐₓ*Δt))
	#εₓ *= nb/(2*ri) .*(sin.(π*fₘₐₓ*ri*Δz/c))./(sin.(π*fₘₐₓ*Δt))
	meh = nb./(2*UR)
	mjj = (2 .-Γ.*Δt)./(2 .+Γ.*Δt)
	mqj = (2 .*Ω.^2 .*Δt)./(2 .+Γ.*Δt)
	mej = (2*Ω.*sqrt.((ε₀.-εₓ)).*Δt)./(2 .+Γ.*Δt)
	mhd = nb/2
	mde = 1./(εₓ)
	mqe = (Ω .*sqrt.((ε₀.-εₓ)))./(εₓ)
	#Initialize fields and Boundary Terms
	E, H, J, Q, D             =zeros(ComplexF64,Nᵣ),zeros(ComplexF64,Nᵣ),zeros(ComplexF64,Nᵣ),zeros(ComplexF64,Nᵣ),zeros(ComplexF64,Nᵣ)
	e₁, e₂ = 0,0
	h₁, h₂ = 0,0

	#Initialize Reflected and Transmitted Fields
	Er = zeros(ComplexF64,Steps)
	Et = zeros(ComplexF64,Steps)

	#Setup Fourier Transforms
	Freq = fftshift(fftfreq(Steps,1/Δt))
	Esf = fftshift(fft(Eₛ))

	#Main FDTD Loop
	for T in 1:Steps
		#Record H at Boundary
		h₂ = h₁
		h₁ = H[1]
		#Update H from E
		for nz in 1:Nᵣ-1
			H[nz] = H[nz] + meh[nz]*(E[nz+1] - E[nz])
		H[Nᵣ] = H[Nᵣ] + meh[Nᵣ]*(e₂-E[Nᵣ])
		#H Source
		H[1] = H[1] - meh[1]*Eₛ[T]
		#Update J from Q and E
		for nz in 1:Nᵣ
			J[nz] = mjj[nz]*J[nz] - mqj[nz]*Q[nz] + mej[nz]*E[nz] 
		#Update Q from J
		for nz in 1:Nᵣ
			Q[nz] = Q[nz] + Δt*J[nz]
		#Update D from H
		D[1] = D[1] + mhd*(H[1] - h₂)
		for nz in 2:Nᵣ
			D[nz] = D[nz] + mhd*(H[nz] - H[nz-1])
		#D Source
		D[2] = D[2] - mhd*Hₛ[T]
		#Record E at Boundary
		e₂ = e₁
		e₁ = E[Nᵣ]
		#Update E from D and Q
		for nz in 1:Nᵣ
			E[nz] = mde[nz]*D[nz]-mqe[nz]*Q[nz]
		#Update reflected and transmitted fields
		Er[T] = E[1]
		Et[T] = E[Nᵣ]

	#Determine optical properties from FDTD values
	Erf = fftshift(fft(Er))
	Etf = fftshift(fft(Et))

	r = Erf./Esf
	tr = Etf./Esf

	Ref = abs.(r).^2
	Trn = abs.(tr).^2
	Con = Ref.+Trn

	n = (1 .-r)./(1 .+r)
	ε₁ = real(n).^2-imag(n).^2
	ε₂ = 2*real(n).*imag(n)

	#Determine optical properties from analytical solution
	εₜₑₛₜ = εₓ[n₁].+(Ω[n₁]^2*(ε₀[n₁]-εₓ[n₁]))./(Ω[n₁]^2 .-Freq.^2 .-im.*Γ[n₁]*Freq) 
	nₜₑₛₜ = sqrt.(εₜₑₛₜ)
	Rₜₑₛₜ=abs.((1 .-nₜₑₛₜ)./(1 .+nₜₑₛₜ)).^2

I would appreciate any help in fixing my code asap as I have not been able to fix this issue for almost a month now.

1 Like