This is likely the physically expected behavior. If your gain is coupling
to a mode of your system whose loss rate is less than that of the rate of
stimulated emission, your system will begin to lase, i.e. the field will
begin to grow exponentially. In real, physical systems, this is then
compensated by gain saturation, which prohibits the fields from growing
exponentially forever. However, as the Lorentzian susceptibility model does
not contain this physics, nothing prohibits the fields from continuing
their exponential growth.

On Sat, Jul 25, 2020 at 2:02 PM 裴延波 <peiya...@163.com> wrote:

> I am trying to use dispersive complex epsilon to describe gain in my
> calculation. The parameters for the epsilon is defined as follows.
> freq_32 = 2            # emission frequency  (units of 2\pi c/a)
> gamma_32 = 0.306    # FWHM emission linewidth in sec^-1 (units of 2\pi c/a)
> sigma_32 = 1e-4
> susceptibilities = [mp.LorentzianSusceptibility(frequency=freq_32,
> gamma=-gamma_32, sigma=sigma_32)]
> geometry = [mp.Block(center=mp.Vector3(z=0),
>                      size=mp.Vector3(mp.inf,mp.inf,dcell),
>
> material=mp.Medium(epsilon=2.25,E_susceptibilities=susceptibilities))]
>
> Here gamma is negative indicating the material has gain, just as that in
> the tutorial. However, when running the field increase to infinity even
> though sigma is very very small. I try to set gamma positive and sigma
> negative. In this case there is a result which looks normal. I don't  know
> why and what is the problem in my code. Can anyone explain this. Thanks in
> advance.
>
> The following is my full python code.
>
> import meep as mp
> import math
> resolution = 100
> dimensions = 3
> ns = 1.0
> nlead = ns
> dlead = 2.0
> npad = ns
> dpad = 2.0
> dpml = 2.0
> Ncell = 20
> dcell = 96
> sz = dcell + dlead + dpad + 2*dpml
> cell_size = mp.Vector3(0,0,sz)
> pml_layers = [mp.PML(dpml)]
> freq_32 = 2            # emission frequency  (units of 2\pi c/a)
> gamma_32 = 0.306    # FWHM emission linewidth in sec^-1 (units of 2\pi c/a)
> df1 = gamma_32/2/math.pi
> sigma_32 = 1e-4      # dipole coupling strength (hbar = 1)
> default_material = mp.Medium(index=ns)
> sources = [mp.Source(mp.GaussianSource(freq_32, fwidth=df1),
> component=mp.Ex, center=mp.Vector3(-dcell/2-dpad/2))]
>
> sim = mp.Simulation(cell_size=cell_size,
>                     sources=sources,
>                     resolution=resolution,
>                     boundary_layers=pml_layers,
>                     dimensions = dimensions,
>                     default_material=default_material)
>
> nfreq = 50 #number of frequencies at which to compute flux
> pt = mp.Vector3(0,0,dcell/2+dpad/2)
> flux_detection_point = mp.FluxRegion(center=pt)
> incidence = sim.add_flux(freq_32,df1,nfreq,flux_detection_point)
> sim.run(until_after_sources=mp.stop_when_fields_decayed(50,mp.Ex,pt,1e-3))
> incident_flux = mp.get_fluxes(incidence)
> sim.reset_meep()
> susceptibilities = [mp.LorentzianSusceptibility(frequency=freq_32,
> gamma=-gamma_32, sigma=sigma_32)]
> geometry = [mp.Block(center=mp.Vector3(z=0),
>                      size=mp.Vector3(mp.inf,mp.inf,dcell),
>
> material=mp.Medium(epsilon=2.25,E_susceptibilities=susceptibilities))]
> geometry.append(mp.Block(center=mp.Vector3(z=-sz/2+(dpml+dlead)/2),
>                          size=mp.Vector3(mp.inf,mp.inf,dpml+dlead),
>                          material=mp.Medium(index=nlead)))
> geometry.append(mp.Block(center=mp.Vector3(z=sz/2-(dpml+dpad)/2),
>                          size=mp.Vector3(mp.inf,mp.inf,dpml+dpad),
>                          material=mp.Medium(index=npad)))
> sim = mp.Simulation(cell_size=cell_size,
>                     sources=sources,
>                     resolution=resolution,
>                     boundary_layers=pml_layers,
>                     geometry=geometry,
>                     dimensions = dimensions,
>                     default_material=default_material)
>
> transmission = sim.add_flux(freq_32,df1,nfreq,flux_detection_point)
> sim.run(until_after_sources=mp.stop_when_fields_decayed(50,mp.Ex,pt,1e-3))
> transmitted_flux = mp.get_fluxes(transmission)
> flux_freqs = mp.get_flux_freqs(transmission)
> data1 = open("lasing.dat",'w')
> for ii in range(0,nfreq):
>      data1.write("%f   %f   %f    %f\n"
> %(flux_freqs[ii],incident_flux[ii],transmitted_flux[ii],transmitted_flux[ii]/incident_flux[ii]))
> data1.close()
>
>
>
>
> _______________________________________________
> meep-discuss mailing list
> meep-discuss@ab-initio.mit.edu
> http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/meep-discuss
_______________________________________________
meep-discuss mailing list
meep-discuss@ab-initio.mit.edu
http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/meep-discuss

Reply via email to