Just want to thank you Dr. Rumpf for sharing your knowledge with the world. I have been able to run a 1D FDTD simulation and it was awesome! I know I still have a long way to go and much to learn about FDTD, but I have to admit I never thought one day I could apply FDTD to simulate EM systems. Thank you so much for your guidance and sharing your expertise with us all. You are a great example for all electrical engineering teachers to follow. Blessings from the Republic of Panama!
Hello Professor I am trying to model a dispersive dielectric material using FDTD to determine its frequency-dependent relative permittivity but I can't get it to work. I tried to modify my working model for the dielectric slab you talked about in the previous lectures (8 and 9). Since now I don't know the refractive index of the material (in fact that is what I want to determine) I took n_max as 1. I should also note that the set of equations I am using is a bit more complicated, I am using the phonon coordinate instead of the polarization but it just replaces the constants in the equations so it shouldn't affect the result I think? I am also still updating the values of H and E at the boundaries which seem to be missing in the loop you gave in this lecture but since the equations for D and H do depend on the values of H and E at the boundaries respectively I think keeping them is justified. Anyway, I wanted to ask you what I should take for n_max (I know that near the resonance the refractive index jumps up to 16+ 16i) in this case and what other things I would have to change to get an accurate result. Thank you in advance for your answer.
Sorry, but I am a little confused. You wrote “trying to model a dispersive dielectric material using FDTD to determine its frequency-dependent relative permittivity.” You will need to know the frequency-dependent permittivity in order to run a simulation. I think maybe you are tweaking the parameters until your simulation matches experiments in order to determine the properties? Setting nmax to 1 when you don’t know the properties is also confusing. Again, you cannot run a simulation without entering some kind of values. In this case, you should know nmax. I do not know what “phonon coordinate” is to help you. Which block diagram are you looking at? It should include recording E and H at the boundaries. From your last statement, choose nmax = 16.
@@empossible1577 Hello once again thank you for your reply. To clarify some things I am using the Lorentz model on dipoles to determine the propagation of light. Here the polarization is simply given by: P=NZQ/sqrt(mu*N) where N is the number of unit cells per unit volume, mu is the effective mass, Z is the charge of the ions of a dipole, and finally, Q is the phonon coordinate given by Q=sqrt(N*mu)*W where W is the effective displacement between the two ions of the dipole. In other words, here a phonon is the different vibrational normal modes of a charged harmonic oscillator (which is the dipole). We can rewrite the governing equation as: partial of J + gamma*J+omega^2*Q=sqrt(N/mu)*Z*E, where J= partial of Q instead of P Taking a monochromatic EM wave (of angular frequency f) the solution for Q is also a monochromatic wave giving: Q(t)=(sqrt(N/mu)*Z)/(omega^2-f^2-i*gamma*f)*E(t) plugging this equation back in for the equation for the polarization we get a similar result as the one you got at 47:48 in the lecture (or slide 24 in the lecture notes). Now using the relation D=epsilon_0*E+P we get: epsilon_r(omega)= 1+(NZ^2)/(epsilon_0*mu)*1/(omega^2-f2-i*gamma*f) =1+(omega^2(epsilon_0-1)/(omega^2-f2-i*gamma*f) which is the same relation you gave at 38:25 (or slide 41 from the lecture notes) once we define omega_p^2 as the coefficient at the start. Now this is the analytical solution for a monochromatic wave. Now using the loop you provided at the end (or slide 52 in the lecture notes which is the one that is missing the steps of updating the values at the boundary) and replacing P with Q (for now this makes no difference except for the coefficient for P which is now Q) in the equation of E we should be able to retrieve this result through FDTD but right now I am having difficulties with this step. That's why I wanted to ask you about the parameters and implementation of this system of equations. I know analytically that the refractive index goes up to 16+16i (but didn't know how to choose nmax to match that). I will try using 16 as per your suggestion. I also wanted to ask you what would happen if I tried to use a higher value just out of curiosity. As for your question on how I can run the simulation without the frequency-dependent relative permittivity, I think epsilon_0 and epsilon_infinity (as the value fo the relative permittivity at 0 and infinite frequency respectively) should be enough the run the simualtion. Like I said I was confused as to how to apply this to the algorithm from before. I should also add that the source I am using is a modulated Gaussian pulse with bandwidth extending from 20 to 32 THz since omega in my case is about 25 THZ. So in this case I am also taking fmax=5 which is not the maximum frequency anymore but just the difference between 20 and 26 THz (which is the central frequency of the modulated Gaussian). So I suppose I would like to ask whether this would pose any problems as well or i should I just change tau to be 0.1 and keep fmax at 32 THz or a slightly higher value. I am using this range as this is the range is slightly bigger than the so-called reststrahlen band where the reflectivity is high. To give my question more context, I plan to add some nonlinear terms to Q afterwards to replicate the results of this paper: www.pnas.org/doi/full/10.1073/pnas.1809725115. However, at this stage, I am excluding the nonlinear terms to make sure I have a working model. Once again thank you for your answer and assistance in advance
@@empossible1577 Hello I wanted to give you an update on my previous comment. After your comment, I tried to put fmax=32 THz and nmax=16 and put tau=0.1 for the modulated Gaussian pulse. Unfortunately, I still can't get accurate results. In fact, weirdly I am retrieving the previous sinusoidal pattern instead of the curve with the restrahlenband (the area where the reflectivity jumps up and then down). So I am asking for your help once again. I have one more question you previously said the spacer region should be as big as the minimum wavelength for 2D and 3D FDTD but is not required for 1D FDTD. Does this still apply to dispersive materials? Thank you for your help and sorry for bothering you once again hopefully I am not asking too much.
@@emredurmus2506 No need for spacer regions in this specific implementation of 1D FDTD regardless of material properties. The reason for the spacer region is that many absorbing boundaries do not handle evanescent fields very well. Inside the perfectly matched layer (PML) absorbing boundary, evanescent waves become propagating waves that become an escape path for power in the simulation that you cannot account for. The spacer region is meant to allow the evanescent fields to sufficient decay before touching the PML. Many times you will hear that the spacer regions need to be a quarter or half wavelength, but this is only a rule-of-thumb. Some devices have very large evanescent fields and these simulation will require much larger spacer regions. This specific implementation of 1D FDTD uses a perfect absorber so it is never a problem. It could can handle evanescent fields without problem. Hope this helps!
@@empossible1577 Yes thank you for your response I just had one more question. For my source, I am using a modulated Gaussian with a central frequency of 26 THz how would I modify fmax in the algorithm in this case I have set tau to 0.1 since I want the range of the source to be between 20 and 30 THz.
Hi Professor, P = e0*xe*E (from slide26, Lorentz-Drude model, 21 century Electromagnetics). Why e0 is missing in the last question of slide 43 at 46:58? Thanks.
Dr. Rumpf, Thank you for your amazing videos. Whenever I am stuck with my research, the first thing I do is to refer to your videos and a couple of books. Recently, we have discovered an error in our research and couldn't find a solution to it. I would appreciate having your views on it. Can I email you about my research and discuss the problem?
Hello, great lectures in explaining the FDTD concept, I am trying to model a three-layer substrate using 1-D simulation and exciting it with a triangular Linear Frequency Continuous Wave modulated sine wave. At the internal boundaries of the substrate where the dielectric permittivity, conductivity and permeability change, should I incorporate empty cells for TF/SF at each boundary or just place the whole substrate in the middle and continue as described in the lecture videos? I only want to look at the reflections that are caused by the incident wave, not the waves which are transmitted through the substrate by internal reflections. Thanks!
The TF/SF is for the source and should not affect how your construct the device on your grid. For simulating reflection, simply fill the first part of your grid with air, then a few more cells for your first layer, then a few more cells for your second, layer, and so on until you can plot the permittivity and see a picture of your three-layer substrate. The TF/SF interface will go in the air region above your three layers.
Hi, I'm trying to incorporate losses into the dielectric that i'm trying to model. For practice, I've modeled a slab of arbitrary conductivity and thickness. For some reason for any value of sigma other than 0 that I give, the slab acts as a perfect conductor (reflects the transmitted wave completely) and I'm not sure why. Any help would be greatly appreciated. Here's the code. ER = [ones(1,100) 4*ones(1,50) ones(1,50)]; UR = ones(1,Nz); sigma = [zeros(1,100) 3*ones(1,50) zeros(1,50)]; const_E1 = (2.*E0.*ER-dt.*sigma)./(2.*E0.*ER+dt.*sigma) const_E2 = (2.*E0.*C0.*dt)./(2.*E0.*ER+dt.*sigma) const_H = (C0*dt)./UR; const_E = (C0*dt)./ER; Ey = zeros(1,Nz); Hx = zeros(1,Nz); e3=0;e2=0;e1=0;h3=0;h2=0;h1=0; NFREQ=1000; FREQ = linspace(0,fmax,NFREQ); k = exp(1i*2*pi*dt*FREQ); REF = zeros(1,NFREQ); TRN = zeros(1,NFREQ); SRC = zeros(1,NFREQ); for t = 1:steps
for nz = 1:Nz-1 Hx(nz) = Hx(nz) + const_H(nz)*( Ey(nz+1)-Ey(nz))/dz; end Hx(Nz) = Hx(Nz) + const_H(Nz)*(e3-Ey(Nz))/dz; Hx(49)=Hx(49)-const_H(49)*E_source(t)/dz; h3=h2; h2=h1; h1=Hx(1);
No. If there are any atomic resonances at a frequency much higher than you are interested in for your simulation, there is no need to model that resonance. Instead, it is only necessary to include the offset due to that resonance. So eps(infinity) is the sum of all the offsets from all the resonances at much higher frequencies. Think of it as a background permittivity.
Dear Prof., For incorporating losses into the simulation, the conductivity is related to the complex permittivity (and index), is that right? So far, the equation that I have found relating both is sigma = omega*eps2/(4*pi) where sigma is the conductivity, omega is the angular frequency, and eps2 is the imaginary part of the complex permittivity. Is this correct? What other ways would you recommend I get the conductivity value from? (Theoretical models or experimental data perhaps?) Thanks!
I usually simulate devices in stages of increasing complexity. Consider first, just setting the electric field to zero where metals exist. This is the PEC approximation. One easy way to do this is to set the permittivity something crazy large like 1e6. The second stage of complexity is incorporating a constant conductivity, usually this is the DC conductivity. There are lots of tables on the internet giving DC conductivity for many different materials. The third stage of complexity is incorporating something a Drude or Lorentz-Drude model to account for the dispersion of conductivity. This is particularly important at very high frequencies. Finding these values is much more challenging. Often, people will measure the conductivity themselves and then do a curve fit to to extract the parameters. The best thing is to search the literature for others that have simulated metals with the Drude model, or whatever model you are using.
Hi, I have 2 more questions, how can I take into account the full anisotropic material with the frequency-dependent model?. should I create nine Lorenzt models for the polarization like simulating each component of the permitivity tensor and then used the interpolation for the off-diagonal componentes (so everything is in the middle point of the curl as exposed in another of your lectures)?. and the second question is; how can I add a fixed magnetic field to the model?, I am trying to get the magneto-optical response when applying a fixed magnetic field in a layered media, so in the analytical solution it will be easy by only adding some fixed components to the off diagonal terms in the permitivity tensor, but in this case I don't know if physically it will be the same by adding these magnetization in the relative pemitivity of the free space (the constant Epsilon infinity showed in the slides). Again Thanks a lot for your time.
First, I suggest you identify the specific physics want to incorporate. For simple anisotropy, there is no need for Lorentz models. The lorentz models enter the picture when you want to account for dispersion. If you were to combine dispersion and anisotropy, you should only use three Lorentz models because we live in a world restricted to three dimensions. You would iterate those models to get the three principle permitivity values. From there, you would rotate the diagonal tensor to calculate the off-diagonal components. In FDTD, there is really no need to incorporate a static field value. That field is modifying the electrical properties of some material. In mnost cases, you can just use the modified material properties and not worry about the static field. If you are looking at layered-media with anisotropy, dispersion, etc., take a look at the transfer matrix method. I think it may serve you better.
@@empossible1577, yes, the transfer matrix works for most of the cases, but now, I am studying layers that have inside nano disk or other types of geometrical "anisotropies" (maybe in that case it is not called layered media). Actually what I need to incorporate is absorbance (and also measure the Kerr or Farady angle but I haven´t get to that yet), but maybe I misunderstood something with the absorption model presented in the lecture witch I believe only allows one frequency at a time, And I though that the next section (Frequecny dependent material), was for implicitly incorporate the absorption with the Lorentz model and in that way stimulate the system with a Gaussian pulse. So if the Lorentz model does not include abortion, I will have another question, If absorption is only manageable with the current therm of the maxwell Equation, I am a little be confused with what will be the value of the sigma therm (the imaginary magnitud of the refractive index of the material?). or is there anyway to have a frequently dependent material that also includes absorption. Again thanks a LOT for your time and patience
@@okar202 The Lorentz model does include absorption. Loss is introduced through the gamma term. It is also possible to introduce simple loss through the conductivity term. I have even seen people do both! If your device size is small, I recommend using a frequency-domain method. It will make incorporating all your crazy stuff very simple. Time-domain will only offer you benefits for nonlinear material properties.
Hi, thanks a lot for the lectures, they are helping me a lot with my thesis, I have a question, when adding the correction for numerical dispersion, the formula involves the wave vector k (that dependes on the frequency) but how can I add the compensation when the system is perpetuated with a Gaussian pulse that has a lot of frequencies, I'm sure I'm missing something but I don't get it. Thanks in advance
Unfortunately, this simple technique can only cancel the dispersion for one frequency in one direction. When you have a spread of frequencies or a spread of angles, all you can do is base the compensation on an average. If the dispersion is a serious problem for you, there are more advanced ways of handling it that give much better results. For examples, you could adopt a hexagonal grid (Lecture 18) or use the M24 algorithm (Lecture 25).
I didn't implement frequency dependent materials. With that said, what should be expected with the reflection of a wave as I increase its frequency? As we know from skin depth equation, a high frequency wave will be attenuated faster than a low frequency wave. However, in my simulations it appears an increase in frequency is doing the other way around, making the wave able to go through the lossy material. Do you know what could be wrong in my equations, or maybe without frequency dependent materials equations the changes in frequency won't work?
The PML is a strange thing as it introduces strong anisotropy, dispersion, and loss. Some dependence on frequency, polarization, and angle of incidence is expected, but should still perform well for you. Without seeing how much variation you are seeing, it is hard to see if what you are getting is correct or not.
using very low conductivity (sigma = 3.54e-3), a 10 MHz wave gets almost completely reflected when reaching the first interface of the material and attenuates the rest before reaching the other end. When I increase the frequency to 10 GHz and further, the wave simply almost ignores the materials, and gets transmitted right through it.
Are you simulating reflection from a PML or are you simulating a structure with conductivity? I do not have enough information to tell what your problem is. Conductivity alters impedance so your wave may be reflected instead of absorbed. Conductivity alters absorption coefficient so if the wave is not reflected it is absorbed. If you are simulating a slab of material, you have resonances that could also be altering the response you are seeing. I suggest simulating something simple like a slab of material that is a quarter-wavelength thick at your center frequency and has some conductivity that is low enough to allow some transmission but high enough to cause some attenuation. This is something easily replicated by another numerical technique like the transfer matrix method.
Respected Sir, What if the dielectric slab has varying refractive index with respect to frequency? then should I calculate the refractive index based constants for every frequency?
+vinpremkumar There are two ways to do this. The least efficient way to do this is to simulate one frequency at a time with a pure frequency source and manually change the refractive index each time. A more efficient way to do this is to incorporate the dispersion directly into your FDTD algorithm. There are many ways to do this, but Z transforms and the Lorentz-Drude seem to be the most popular. See the last section in Lecture 10.
Hi, I am trying to include dispersion in the FDTD model. So, What should be the intial value for P(at time t) and J(at time t-dt) while including dispersion in the FDTD model. Assuming there is only Py component Because as per the lecture for future values of J(t+dt), we need J(t-dt) ,P(t) and E(t).
Hi , I am facing issue while implementing dispersion in the 1D FDTD Below is the additional/modified line of code I have added in the original code to include dispersion. (Made some correction but still code is unstable) %SLAB PROPERTIES dslab=0.0762; erair=1; erslab=12; Gamma=1.2e14; omega_t=9.65e14; epsilon_dc=12; param=e0*(epsilon_dc-1)*omega_t^2; %Here I am assuming fm=1 %-----other code------ % INITIALIZE MATERIALS TO FREE SPACE ER=f*ones(1,Nz); UR=f*ones(1,Nz); GAMMA_0=zeros(1,Nz); OMEGA_t=zeros(1,Nz); PARAM_t=zeros(1,Nz); %----other code------- %ADD SLAB ER(nz1:nz2)=f*erslab; GAMMA_0(nz1:nz2)=Gamma; OMEGA_t(nz1:nz2)=omega_t; PARAM_t(nz1:nz2)=param; C1=(2-GAMMA_0*dt)./(2+GAMMA_0*dt); C2=2*OMEGA_t.^2.*dt./(2+GAMMA_0*dt); C3=2*dt*PARAM_t./(2+GAMMA_0*dt); % assuming fm=1 %-------other code------ %INITILIZE FIELDS Ey=zeros(1,Nz); Hx=zeros(1,Nz); Py=zeros(1,Nz); Jy=zeros(1,Nz); Dy=zeros(1,Nz); %-------other code------ %IN THE MAIN LOOP for nz=1:Nz-1 Hx(nz)=Hx(nz)+mHx(nz)*(Ey(nz+1)-Ey(nz))/dz; end Hx(Nz)=Hx(Nz)+mHx(Nz)*( E3 -Ey(Nz))/dz; % H source Hx(nz_src-1) = Hx(nz_src-1)-mHx(nz_src-1)*Esrc(T)/dz; % Record H-field at boundary H3=H2; H2=H1; H1=Hx(1);
%Update J from P and E for nz=1:Nz Jy(nz)=C1(nz)*Jy(nz)-C2(nz)*Py(nz)+C3(nz)*Ey(nz); end
%Update P from J for nz=1:Nz Py(nz)=Py(nz)+dt*Jy(nz); end
% Update D from H Dy(1)=Dy(1)+e0*c0*dt*(Hx(1)- H3 )/dz; for nz =2 : Nz Dy(nz)=Dy(nz)+e0*c0*dt*(Hx(nz)-Hx(nz-1) )/dz; end % D source Dy(nz_src)=Dy(nz_src)-e0*c0*dt*Hsrc(T)/dz; %Update E from D for nz=1:Nz Ey(nz)=1./e0*(Dy(nz)-Py(nz)); end But still the code is very unstable. The output is unbounded. :( I will be very thankful for any kind of advice/suggestion in this regard Once again thanks for such an insightful and awesome lectures. Usually,I find it very difficult to learn FDTD from the book. but your lectures really helps me to understand the basics of FDTD as well as how to implement in the simulation.
I hope you are doing great Dr. Rumpf. I would like to know what is the problem in the code above too as I am facing the same problem. The quantities become unbounded very quickly. I have the polarization terms derived in a different way thought. And it works. By 'works' I mean it doesn't produce unbounded outputs. But there is no way to tell if the code is actually right. Conservation of energy method of checking will not work here because there is absorption here. If there were some practical examples like lecture 9 we could have verified our results.
Hello !! I made a disperssion code with matlab.Does the fact that at some points conservation (sum of R and T) is below 1 derives from the material's disperssion loss?
If the dispersion you incorporated includes loss, this is certainly expected. When I am simulating things with loss, I will sometimes turn off of the loss, run the simulation, and look for conservation to be a little more sure things are correct. Then I turn the loss back on and run the simulation.
I have made a code for 1FDTD with cole-cole media .I have an analytical formula for the Reflectance so i can see if i made the code right.The fact is that the Fourier of FDTD is the same with the analytical results only when the material is only a few cells .When i put more cells ,the reflectance stays at a price tha doesnt change if i add more, but its different from analytical for 10% .Of course it changes with frequency always ,my problem is with the length of the material.Is it acceptable or do u believe i made a mistake ? The analytical formula is the : abs[R(w)]=abs [ (1-sqrt(epsilon_r(w))) / (1-sqrt(epsilon_r(w))) ]
Are you simulating transmission through a slab of material? The equation you have does not have thickness as a parameter so I am not sure what exactly you are modeling. Can you clarify?
Make sure you are always simulating a slab of the same physical size. Multiple the number of cells wide by the grid resolution parameter dz to check the thickness. Be sure that is consistent.
Just want to thank you Dr. Rumpf for sharing your knowledge with the world. I have been able to run a 1D FDTD simulation and it was awesome! I know I still have a long way to go and much to learn about FDTD, but I have to admit I never thought one day I could apply FDTD to simulate EM systems. Thank you so much for your guidance and sharing your expertise with us all. You are a great example for all electrical engineering teachers to follow. Blessings from the Republic of Panama!
Greeting Republic of Panama!! This is great to hear! Congrats on getting your code to work!
Hello Professor I am trying to model a dispersive dielectric material using FDTD to determine its frequency-dependent relative permittivity but I can't get it to work. I tried to modify my working model for the dielectric slab you talked about in the previous lectures (8 and 9).
Since now I don't know the refractive index of the material (in fact that is what I want to determine) I took n_max as 1. I should also note that the set of equations I am using is a bit more complicated, I am using the phonon coordinate instead of the polarization but it just replaces the constants in the equations so it shouldn't affect the result I think? I am also still updating the values of H and E at the boundaries which seem to be missing in the loop you gave in this lecture but since the equations for D and H do depend on the values of H and E at the boundaries respectively I think keeping them is justified.
Anyway, I wanted to ask you what I should take for n_max (I know that near the resonance the refractive index jumps up to 16+ 16i) in this case and what other things I would have to change to get an accurate result.
Thank you in advance for your answer.
Sorry, but I am a little confused. You wrote “trying to model a dispersive dielectric material using FDTD to determine its frequency-dependent relative permittivity.” You will need to know the frequency-dependent permittivity in order to run a simulation. I think maybe you are tweaking the parameters until your simulation matches experiments in order to determine the properties?
Setting nmax to 1 when you don’t know the properties is also confusing. Again, you cannot run a simulation without entering some kind of values. In this case, you should know nmax.
I do not know what “phonon coordinate” is to help you.
Which block diagram are you looking at? It should include recording E and H at the boundaries.
From your last statement, choose nmax = 16.
@@empossible1577 Hello once again thank you for your reply. To clarify some things I am using the Lorentz model on dipoles to determine the propagation of light. Here the polarization is simply given by:
P=NZQ/sqrt(mu*N) where N is the number of unit cells per unit volume, mu is the effective mass, Z is the charge of the ions of a dipole, and finally, Q is the phonon coordinate given by Q=sqrt(N*mu)*W where W is the effective displacement between the two ions of the dipole. In other words, here a phonon is the different vibrational normal modes of a charged harmonic oscillator (which is the dipole). We can rewrite the governing equation as:
partial of J + gamma*J+omega^2*Q=sqrt(N/mu)*Z*E, where J= partial of Q instead of P
Taking a monochromatic EM wave (of angular frequency f) the solution for Q is also a monochromatic wave giving:
Q(t)=(sqrt(N/mu)*Z)/(omega^2-f^2-i*gamma*f)*E(t)
plugging this equation back in for the equation for the polarization we get a similar result as the one you got at 47:48 in the lecture (or slide 24 in the lecture notes). Now using the relation D=epsilon_0*E+P we get:
epsilon_r(omega)= 1+(NZ^2)/(epsilon_0*mu)*1/(omega^2-f2-i*gamma*f)
=1+(omega^2(epsilon_0-1)/(omega^2-f2-i*gamma*f)
which is the same relation you gave at 38:25 (or slide 41 from the lecture notes) once we define omega_p^2 as the coefficient at the start. Now this is the analytical solution for a monochromatic wave. Now using the loop you provided at the end (or slide 52 in the lecture notes which is the one that is missing the steps of updating the values at the boundary) and replacing P with Q (for now this makes no difference except for the coefficient for P which is now Q) in the equation of E we should be able to retrieve this result through FDTD but right now I am having difficulties with this step. That's why I wanted to ask you about the parameters and implementation of this system of equations. I know analytically that the refractive index goes up to 16+16i (but didn't know how to choose nmax to match that). I will try using 16 as per your suggestion. I also wanted to ask you what would happen if I tried to use a higher value just out of curiosity.
As for your question on how I can run the simulation without the frequency-dependent relative permittivity, I think epsilon_0 and epsilon_infinity (as the value fo the relative permittivity at 0 and infinite frequency respectively) should be enough the run the simualtion. Like I said I was confused as to how to apply this to the algorithm from before. I should also add that the source I am using is a modulated Gaussian pulse with bandwidth extending from 20 to 32 THz since omega in my case is about 25 THZ. So in this case I am also taking fmax=5 which is not the maximum frequency anymore but just the difference between 20 and 26 THz (which is the central frequency of the modulated Gaussian). So I suppose I would like to ask whether this would pose any problems as well or i should I just change tau to be 0.1 and keep fmax at 32 THz or a slightly higher value. I am using this range as this is the range is slightly bigger than the so-called reststrahlen band where the reflectivity is high.
To give my question more context, I plan to add some nonlinear terms to Q afterwards to replicate the results of this paper: www.pnas.org/doi/full/10.1073/pnas.1809725115. However, at this stage, I am excluding the nonlinear terms to make sure I have a working model.
Once again thank you for your answer and assistance in advance
@@empossible1577 Hello I wanted to give you an update on my previous comment. After your comment, I tried to put fmax=32 THz and nmax=16 and put tau=0.1 for the modulated Gaussian pulse. Unfortunately, I still can't get accurate results. In fact, weirdly I am retrieving the previous sinusoidal pattern instead of the curve with the restrahlenband (the area where the reflectivity jumps up and then down). So I am asking for your help once again. I have one more question you previously said the spacer region should be as big as the minimum wavelength for 2D and 3D FDTD but is not required for 1D FDTD. Does this still apply to dispersive materials?
Thank you for your help and sorry for bothering you once again hopefully I am not asking too much.
@@emredurmus2506 No need for spacer regions in this specific implementation of 1D FDTD regardless of material properties.
The reason for the spacer region is that many absorbing boundaries do not handle evanescent fields very well. Inside the perfectly matched layer (PML) absorbing boundary, evanescent waves become propagating waves that become an escape path for power in the simulation that you cannot account for. The spacer region is meant to allow the evanescent fields to sufficient decay before touching the PML. Many times you will hear that the spacer regions need to be a quarter or half wavelength, but this is only a rule-of-thumb. Some devices have very large evanescent fields and these simulation will require much larger spacer regions.
This specific implementation of 1D FDTD uses a perfect absorber so it is never a problem. It could can handle evanescent fields without problem.
Hope this helps!
@@empossible1577 Yes thank you for your response I just had one more question. For my source, I am using a modulated Gaussian with a central frequency of 26 THz how would I modify fmax in the algorithm in this case I have set tau to 0.1 since I want the range of the source to be between 20 and 30 THz.
Hi Professor, P = e0*xe*E (from slide26, Lorentz-Drude model, 21 century Electromagnetics). Why e0 is missing in the last question of slide 43 at 46:58? Thanks.
It was normalized out.
@@empossible1577 Thanks for your reply. But if it was normalized out, why e0 appeared in free space response term of D=e0*e(infinity)*E + P at 46:58?
HI Professor, at 52:07, why it's handling D TF/SF source, instead of E TF/SF source? Is it ok to handle E TF/SF source as previous? Thanks.
I incorporate the source from the curl operation. The curl operation happens when calculating D, thus that is where I put the TF/SF correction step.
Dr. Rumpf,
Thank you for your amazing videos. Whenever I am stuck with my research, the first thing I do is to refer to your videos and a couple of books. Recently, we have discovered an error in our research and couldn't find a solution to it. I would appreciate having your views on it. Can I email you about my research and discuss the problem?
Thank you for the compliment! I am happy the course materials are helping you. You are certainly welcome to e-mail me.
@@empossible1577 Thank you for your response, Dr. Rumpf. I will send you the email in a couple of minutes.
Hello, great lectures in explaining the FDTD concept, I am trying to model a three-layer substrate using 1-D simulation and exciting it with a triangular Linear Frequency Continuous Wave modulated sine wave. At the internal boundaries of the substrate where the dielectric permittivity, conductivity and permeability change, should I incorporate empty cells for TF/SF at each boundary or just place the whole substrate in the middle and continue as described in the lecture videos? I only want to look at the reflections that are caused by the incident wave, not the waves which are transmitted through the substrate by internal reflections.
Thanks!
The TF/SF is for the source and should not affect how your construct the device on your grid. For simulating reflection, simply fill the first part of your grid with air, then a few more cells for your first layer, then a few more cells for your second, layer, and so on until you can plot the permittivity and see a picture of your three-layer substrate. The TF/SF interface will go in the air region above your three layers.
Hi, I'm trying to incorporate losses into the dielectric that i'm trying to model. For practice, I've modeled a slab of arbitrary conductivity and thickness. For some reason for any value of sigma other than 0 that I give, the slab acts as a perfect conductor (reflects the transmitted wave completely) and I'm not sure why. Any help would be greatly appreciated. Here's the code.
ER = [ones(1,100) 4*ones(1,50) ones(1,50)];
UR = ones(1,Nz);
sigma = [zeros(1,100) 3*ones(1,50) zeros(1,50)];
const_E1 = (2.*E0.*ER-dt.*sigma)./(2.*E0.*ER+dt.*sigma)
const_E2 = (2.*E0.*C0.*dt)./(2.*E0.*ER+dt.*sigma)
const_H = (C0*dt)./UR;
const_E = (C0*dt)./ER;
Ey = zeros(1,Nz);
Hx = zeros(1,Nz);
e3=0;e2=0;e1=0;h3=0;h2=0;h1=0;
NFREQ=1000;
FREQ = linspace(0,fmax,NFREQ);
k = exp(1i*2*pi*dt*FREQ);
REF = zeros(1,NFREQ);
TRN = zeros(1,NFREQ);
SRC = zeros(1,NFREQ);
for t = 1:steps
for nz = 1:Nz-1
Hx(nz) = Hx(nz) + const_H(nz)*( Ey(nz+1)-Ey(nz))/dz;
end
Hx(Nz) = Hx(Nz) + const_H(Nz)*(e3-Ey(Nz))/dz;
Hx(49)=Hx(49)-const_H(49)*E_source(t)/dz;
h3=h2; h2=h1; h1=Hx(1);
Ey(1) = const_E1(1)*Ey(1) + const_E2(1)*(Hx(1) - h3)/dz;
for nz = 2:Nz
Ey(nz) = const_E1(nz)*Ey(nz) + const_E2(nz)*(Hx(nz) -Hx(nz-1))/dz;
end
Ey(50) = Ey(50)-const_E(50)*H_source(t)/dz;
e3=e2; e2=e1; e1= Ey(Nz);
draw1d(ER,Ey,Hx,dz);
xlim([dz Nz*dz];
xlabel('z');
title(['Field at step' num2str(t) 'of' num2str(steps)]);
drawnow;
for nf=1:NFREQ
REF(nf)=REF(nf)+(k(nf)^t)*Ey(1);
TRN(nf)=TRN(nf)+(k(nf)^t)*Ey(Nz);
SRC(nf)=SRC(nf)+(k(nf)^t)*E_source(t);
end
end
Hi Professor, at 50:51, the value of epsilon(infinity) in the equation is always 1, right? Thanks.
No. If there are any atomic resonances at a frequency much higher than you are interested in for your simulation, there is no need to model that resonance. Instead, it is only necessary to include the offset due to that resonance. So eps(infinity) is the sum of all the offsets from all the resonances at much higher frequencies. Think of it as a background permittivity.
Dear Prof.,
For incorporating losses into the simulation, the conductivity is related to the complex permittivity (and index), is that right? So far, the equation that I have found relating both is
sigma = omega*eps2/(4*pi)
where sigma is the conductivity, omega is the angular frequency, and eps2 is the imaginary part of the complex permittivity. Is this correct?
What other ways would you recommend I get the conductivity value from? (Theoretical models or experimental data perhaps?)
Thanks!
I usually simulate devices in stages of increasing complexity. Consider first, just setting the electric field to zero where metals exist. This is the PEC approximation. One easy way to do this is to set the permittivity something crazy large like 1e6. The second stage of complexity is incorporating a constant conductivity, usually this is the DC conductivity. There are lots of tables on the internet giving DC conductivity for many different materials. The third stage of complexity is incorporating something a Drude or Lorentz-Drude model to account for the dispersion of conductivity. This is particularly important at very high frequencies. Finding these values is much more challenging. Often, people will measure the conductivity themselves and then do a curve fit to to extract the parameters. The best thing is to search the literature for others that have simulated metals with the Drude model, or whatever model you are using.
Hi, I have 2 more questions, how can I take into account the full anisotropic material with the frequency-dependent model?. should I create nine Lorenzt models for the polarization like simulating each component of the permitivity tensor and then used the interpolation for the off-diagonal componentes (so everything is in the middle point of the curl as exposed in another of your lectures)?. and the second question is; how can I add a fixed magnetic field to the model?, I am trying to get the magneto-optical response when applying a fixed magnetic field in a layered media, so in the analytical solution it will be easy by only adding some fixed components to the off diagonal terms in the permitivity tensor, but in this case I don't know if physically it will be the same by adding these magnetization in the relative pemitivity of the free space (the constant Epsilon infinity showed in the slides). Again Thanks a lot for your time.
First, I suggest you identify the specific physics want to incorporate. For simple anisotropy, there is no need for Lorentz models. The lorentz models enter the picture when you want to account for dispersion.
If you were to combine dispersion and anisotropy, you should only use three Lorentz models because we live in a world restricted to three dimensions. You would iterate those models to get the three principle permitivity values. From there, you would rotate the diagonal tensor to calculate the off-diagonal components.
In FDTD, there is really no need to incorporate a static field value. That field is modifying the electrical properties of some material. In mnost cases, you can just use the modified material properties and not worry about the static field.
If you are looking at layered-media with anisotropy, dispersion, etc., take a look at the transfer matrix method. I think it may serve you better.
@@empossible1577, yes, the transfer matrix works for most of the cases, but now, I am studying layers that have inside nano disk or other types of geometrical "anisotropies" (maybe in that case it is not called layered media). Actually what I need to incorporate is absorbance (and also measure the Kerr or Farady angle but I haven´t get to that yet), but maybe I misunderstood something with the absorption model presented in the lecture witch I believe only allows one frequency at a time, And I though that the next section (Frequecny dependent material), was for implicitly incorporate the absorption with the Lorentz model and in that way stimulate the system with a Gaussian pulse.
So if the Lorentz model does not include abortion, I will have another question, If absorption is only manageable with the current therm of the maxwell Equation, I am a little be confused with what will be the value of the sigma therm (the imaginary magnitud of the refractive index of the material?). or is there anyway to have a frequently dependent material that also includes absorption.
Again thanks a LOT for your time and patience
@@okar202 The Lorentz model does include absorption. Loss is introduced through the gamma term. It is also possible to introduce simple loss through the conductivity term. I have even seen people do both!
If your device size is small, I recommend using a frequency-domain method. It will make incorporating all your crazy stuff very simple. Time-domain will only offer you benefits for nonlinear material properties.
Hi, thanks a lot for the lectures, they are helping me a lot with my thesis, I have a question, when adding the correction for numerical dispersion, the formula involves the wave vector k (that dependes on the frequency) but how can I add the compensation when the system is perpetuated with a Gaussian pulse that has a lot of frequencies, I'm sure I'm missing something but I don't get it. Thanks in advance
Unfortunately, this simple technique can only cancel the dispersion for one frequency in one direction. When you have a spread of frequencies or a spread of angles, all you can do is base the compensation on an average. If the dispersion is a serious problem for you, there are more advanced ways of handling it that give much better results. For examples, you could adopt a hexagonal grid (Lecture 18) or use the M24 algorithm (Lecture 25).
I didn't implement frequency dependent materials. With that said, what should be expected with the reflection of a wave as I increase its frequency?
As we know from skin depth equation, a high frequency wave will be attenuated faster than a low frequency wave. However, in my simulations it appears an increase in frequency is doing the other way around, making the wave able to go through the lossy material. Do you know what could be wrong in my equations, or maybe without frequency dependent materials equations the changes in frequency won't work?
The PML is a strange thing as it introduces strong anisotropy, dispersion, and loss. Some dependence on frequency, polarization, and angle of incidence is expected, but should still perform well for you. Without seeing how much variation you are seeing, it is hard to see if what you are getting is correct or not.
using very low conductivity (sigma = 3.54e-3), a 10 MHz wave gets almost completely reflected when reaching the first interface of the material and attenuates the rest before reaching the other end. When I increase the frequency to 10 GHz and further, the wave simply almost ignores the materials, and gets transmitted right through it.
Are you simulating reflection from a PML or are you simulating a structure with conductivity? I do not have enough information to tell what your problem is. Conductivity alters impedance so your wave may be reflected instead of absorbed. Conductivity alters absorption coefficient so if the wave is not reflected it is absorbed. If you are simulating a slab of material, you have resonances that could also be altering the response you are seeing.
I suggest simulating something simple like a slab of material that is a quarter-wavelength thick at your center frequency and has some conductivity that is low enough to allow some transmission but high enough to cause some attenuation. This is something easily replicated by another numerical technique like the transfer matrix method.
Respected Sir,
What if the dielectric slab has varying refractive index with respect to frequency? then should I calculate the refractive index based constants for every frequency?
+vinpremkumar There are two ways to do this. The least efficient way to do this is to simulate one frequency at a time with a pure frequency source and manually change the refractive index each time. A more efficient way to do this is to incorporate the dispersion directly into your FDTD algorithm. There are many ways to do this, but Z transforms and the Lorentz-Drude seem to be the most popular. See the last section in Lecture 10.
+vinpremkumar Thank you. I hope further questions will not be a botheration.
Hi, I am trying to include dispersion in the FDTD model.
So, What should be the intial value for P(at time t) and J(at time t-dt) while including dispersion in the FDTD model. Assuming there is only Py component
Because as per the lecture for future values of J(t+dt), we need J(t-dt) ,P(t) and E(t).
Set everything to zero at the start of the simulation.
Hi , I am facing issue while implementing dispersion in the 1D FDTD
Below is the additional/modified line of code I have added in the original code to include dispersion.
(Made some correction but still code is unstable)
%SLAB PROPERTIES
dslab=0.0762;
erair=1;
erslab=12;
Gamma=1.2e14;
omega_t=9.65e14;
epsilon_dc=12;
param=e0*(epsilon_dc-1)*omega_t^2;
%Here I am assuming fm=1
%-----other code------
% INITIALIZE MATERIALS TO FREE SPACE
ER=f*ones(1,Nz);
UR=f*ones(1,Nz);
GAMMA_0=zeros(1,Nz);
OMEGA_t=zeros(1,Nz);
PARAM_t=zeros(1,Nz);
%----other code-------
%ADD SLAB
ER(nz1:nz2)=f*erslab;
GAMMA_0(nz1:nz2)=Gamma;
OMEGA_t(nz1:nz2)=omega_t;
PARAM_t(nz1:nz2)=param;
C1=(2-GAMMA_0*dt)./(2+GAMMA_0*dt);
C2=2*OMEGA_t.^2.*dt./(2+GAMMA_0*dt);
C3=2*dt*PARAM_t./(2+GAMMA_0*dt); % assuming fm=1
%-------other code------
%INITILIZE FIELDS
Ey=zeros(1,Nz);
Hx=zeros(1,Nz);
Py=zeros(1,Nz);
Jy=zeros(1,Nz);
Dy=zeros(1,Nz);
%-------other code------
%IN THE MAIN LOOP
for nz=1:Nz-1
Hx(nz)=Hx(nz)+mHx(nz)*(Ey(nz+1)-Ey(nz))/dz;
end
Hx(Nz)=Hx(Nz)+mHx(Nz)*( E3 -Ey(Nz))/dz;
% H source
Hx(nz_src-1) = Hx(nz_src-1)-mHx(nz_src-1)*Esrc(T)/dz;
% Record H-field at boundary
H3=H2;
H2=H1;
H1=Hx(1);
%Update J from P and E
for nz=1:Nz
Jy(nz)=C1(nz)*Jy(nz)-C2(nz)*Py(nz)+C3(nz)*Ey(nz);
end
%Update P from J
for nz=1:Nz
Py(nz)=Py(nz)+dt*Jy(nz);
end
% Update D from H
Dy(1)=Dy(1)+e0*c0*dt*(Hx(1)- H3 )/dz;
for nz =2 : Nz
Dy(nz)=Dy(nz)+e0*c0*dt*(Hx(nz)-Hx(nz-1) )/dz;
end
% D source
Dy(nz_src)=Dy(nz_src)-e0*c0*dt*Hsrc(T)/dz;
%Update E from D
for nz=1:Nz
Ey(nz)=1./e0*(Dy(nz)-Py(nz));
end
But still the code is very unstable. The output is unbounded. :(
I will be very thankful for any kind of advice/suggestion in this regard
Once again thanks for such an insightful and awesome lectures. Usually,I find it very difficult to learn FDTD from the book. but your lectures really helps me to understand the basics of FDTD as well as how to implement in the simulation.
I hope you are doing great Dr. Rumpf.
I would like to know what is the problem in the code above too as I am facing the same problem. The quantities become unbounded very quickly. I have the polarization terms derived in a different way thought. And it works. By 'works' I mean it doesn't produce unbounded outputs. But there is no way to tell if the code is actually right. Conservation of energy method of checking will not work here because there is absorption here. If there were some practical examples like lecture 9 we could have verified our results.
Hello !! I made a disperssion code with matlab.Does the fact that at some points conservation (sum of R and T) is below 1 derives from the material's disperssion loss?
If the dispersion you incorporated includes loss, this is certainly expected. When I am simulating things with loss, I will sometimes turn off of the loss, run the simulation, and look for conservation to be a little more sure things are correct. Then I turn the loss back on and run the simulation.
I have made a code for 1FDTD with cole-cole media .I have an analytical formula for the Reflectance so i can see if i made the code right.The fact is that the Fourier of FDTD is the same with the analytical results only when the material is only a few cells .When i put more cells ,the reflectance stays at a price tha doesnt change if i add more, but its different from analytical for 10% .Of course it changes with frequency always ,my problem is with the length of the material.Is it acceptable or do u believe i made a mistake ? The analytical formula is the :
abs[R(w)]=abs [ (1-sqrt(epsilon_r(w))) / (1-sqrt(epsilon_r(w))) ]
Are you simulating transmission through a slab of material? The equation you have does not have thickness as a parameter so I am not sure what exactly you are modeling. Can you clarify?
CEM Lectures i have only a slab .its exactily like the example u have in matlab but instead of e_r constant i have e_r(w) from cole-cole equation .
Make sure you are always simulating a slab of the same physical size. Multiple the number of cells wide by the grid resolution parameter dz to check the thickness. Be sure that is consistent.