In this example, we consider wave propagation through a simple 1d nonlinear medium with a non-zero Kerr susceptibility $\chi^{(3)}$. See also Materials and Units and Nonlinearity. We send in a narrow-band pulse at a frequency $\omega$, and because of the nonlinearity we also get a signal at a frequency $3\omega$.
Since this is a 1d calculation, we could implement it via a 2d cell of Vector3(S,0,0)
, specifying periodic boundary conditions in the $y$ direction. However, this is slightly inefficient since the $y$ periodic boundaries are implemented internally via extra "ghost pixels" in the $y$ direction. Instead, Meep has special support for 1d simulations in the $z$ direction. To use this, we must explicitly set dimensions
to 1 in the Simulation
object constructor, and in that case we can only use $E_x$ (and $D_x$) and $H_y$ field components. This involves no loss of generality because of the symmetry of the problem.
We are going to fill the entire computational cell with the nonlinear medium, so we don't need to use any objects. We can just use the special default_material
which is ordinarily vacuum. PMLs are placed at the $\pm 1$ boundaries. The source will be a Gaussian pulse of $J_x$ just next to the $-z$ PML layer. Since this is a nonlinear calculation, we may want to play with the amplitude of the current/field, so we set the amplitude
property explicitly to a parameter amp
.
from typing import Tuple, List, Union
import numpy as np
from matplotlib import pyplot as plt
import meep as mp
Using MPI version 3.1, 1 processes
Invalid MIT-MAGIC-COOKIE-1 key
def third_harmonic_generation(
k: float, amp: float = 1.0, nfreq: int = 10, flux_spectrum: bool = True
) -> Union[Tuple[List[float], List[float]], Tuple[float, float]]:
"""Computes the transmission spectrum of a plane wave propagating
through a Kerr medium.
Args:
k: strength of Kerr susceptibility.
amp: amplitude of the incident planewave.
nfreq: number of frequencies in flux spectrum.
flux_spectrum: compute the flux spectrum over broad bandwidth (True) or
just the two harmonic frequencies ω and 3ω (False).
Returns:
The frequencies and transmitted flux over the flux spectrum or
the transmitted flux at the harmonic frequencies ω and 3ω.
"""
sz = 100 # size of cell in z direction
fcen = 1 / 3.0 # center frequency of source
df = fcen / 20.0 # frequency width of source
dpml = 1.0 # PML thickness
dimensions = 1
cell = mp.Vector3(0, 0, sz)
pml_layers = [mp.PML(dpml)]
resolution = 25
default_material = mp.Medium(index=1, chi3=k)
sources = [
mp.Source(
mp.GaussianSource(fcen, fwidth=df),
component=mp.Ex,
center=mp.Vector3(0, 0, -0.5 * sz + dpml),
amplitude=amp,
)
]
# frequency range for flux calculation
fmin = fcen / 2.0
fmax = fcen * 4
sim = mp.Simulation(
cell_size=cell,
sources=sources,
boundary_layers=pml_layers,
default_material=default_material,
resolution=resolution,
dimensions=dimensions,
)
mon_pt = mp.Vector3(0, 0, 0.5 * sz - dpml - 0.5)
if flux_spectrum:
trans = sim.add_flux(
0.5 * (fmin + fmax),
fmax - fmin,
nfreq,
mp.FluxRegion(mon_pt),
)
else:
trans1 = sim.add_flux(fcen, 0, 1, mp.FluxRegion(mon_pt))
trans3 = sim.add_flux(
3 * fcen, 0, 1, mp.FluxRegion(mon_pt)
)
sim.run(until_after_sources=mp.stop_when_fields_decayed(50, mp.Ex, mon_pt, 1e-6))
if flux_spectrum:
freqs = mp.get_flux_freqs(trans)
trans_flux = mp.get_fluxes(trans)
return freqs, trans_flux
else:
print(
f"harmonics:, {k}, {amp}, {mp.get_fluxes(trans1)[0]}, "
f"{mp.get_fluxes(trans3)[0]}"
)
return mp.get_fluxes(trans1)[0], mp.get_fluxes(trans3)[0]
Note DFT decimation is turned off by default whenever nonlinearities are present.
In the first part of this tutorial, we plot the transmitted power spectrum for various values of $\chi^{(3)}$. In a linear problem, we normally look at the spectrum over the same frequency range as our source, because other frequencies are zero. In this case, however, we will look from fcen/2
to 4*fcen
, to be sure that we can see the third-harmonic frequency.
nfreq = 400
logk = range(-3, 1)
tflux = np.zeros((nfreq, len(logk)))
for i, lk in enumerate(logk):
freqs, tflux[:, i] = third_harmonic_generation(
10**lk, nfreq=nfreq, flux_spectrum=True
)
----------- Initializing structure... time for choose_chunkdivision = 0.000160748 s Working in 1D dimensions. Computational cell is 0 x 0 x 100 with resolution 25 time for set_epsilon = 0.00186399 s ----------- field decay(t = 100.04): 5.345069974777698e-12 / 5.345069974777698e-12 = 1.0 field decay(t = 150.06): 1.0169446296713276e-08 / 1.0169446296713276e-08 = 1.0 field decay(t = 200.08): 4.650911893321058e-06 / 4.650911893321058e-06 = 1.0 field decay(t = 250.1): 0.0005451205008237801 / 0.0005451205008237801 = 1.0 field decay(t = 300.12): 0.01777078053576906 / 0.01777078053576906 = 1.0 field decay(t = 350.14): 0.13039636481724018 / 0.13039636481724018 = 1.0 field decay(t = 400.16): 0.24555052973214989 / 0.24555052973214989 = 1.0 field decay(t = 450.18): 0.24496560125506897 / 0.24555052973214989 = 0.9976178895736085 field decay(t = 500.2): 0.11044114777393739 / 0.24555052973214989 = 0.4497695358035196 field decay(t = 550.22): 0.012780456545355666 / 0.24555052973214989 = 0.052048173381245705 field decay(t = 600.24): 0.00037674024813005836 / 0.24555052973214989 = 0.0015342677066956937 field decay(t = 650.26): 2.4030087682351996e-06 / 0.24555052973214989 = 9.786208854268963e-06 field decay(t = 700.28): 4.467052624211472e-09 / 0.24555052973214989 = 1.8191989359926036e-08 run 0 finished at t = 700.28 (35014 timesteps) ----------- Initializing structure... time for choose_chunkdivision = 8.3215e-05 s Working in 1D dimensions. Computational cell is 0 x 0 x 100 with resolution 25 time for set_epsilon = 0.000742635 s ----------- field decay(t = 100.04): 5.345069974770661e-12 / 5.345069974770661e-12 = 1.0 field decay(t = 150.06): 1.0169446225969908e-08 / 1.0169446225969908e-08 = 1.0 field decay(t = 200.08): 4.650898281484011e-06 / 4.650898281484011e-06 = 1.0 field decay(t = 250.1): 0.000544949642624633 / 0.000544949642624633 = 1.0 field decay(t = 300.12): 0.017569355547691595 / 0.017569355547691595 = 1.0 field decay(t = 350.14): 0.12281472609090456 / 0.12281472609090456 = 1.0 field decay(t = 400.16): 0.22531439709826295 / 0.22531439709826295 = 1.0 field decay(t = 450.18): 0.22482927464433763 / 0.22531439709826295 = 0.9978469087631637 field decay(t = 500.2): 0.1079829572060197 / 0.22531439709826295 = 0.47925458202711624 field decay(t = 550.22): 0.01266950106036936 / 0.22531439709826295 = 0.05623032182379363 field decay(t = 600.24): 0.00037665107710744167 / 0.22531439709826295 = 0.0016716689299848804 field decay(t = 650.26): 2.403005349498868e-06 / 0.22531439709826295 = 1.0665121183760316e-05 field decay(t = 700.28): 4.467048459368641e-09 / 0.22531439709826295 = 1.982584564900438e-08 run 0 finished at t = 700.28 (35014 timesteps) ----------- Initializing structure... time for choose_chunkdivision = 8.3803e-05 s Working in 1D dimensions. Computational cell is 0 x 0 x 100 with resolution 25 time for set_epsilon = 0.000510164 s ----------- field decay(t = 100.04): 5.345069974700306e-12 / 5.345069974700306e-12 = 1.0 field decay(t = 150.06): 1.0169445518535415e-08 / 1.0169445518535415e-08 = 1.0 field decay(t = 200.08): 4.650762120863413e-06 / 4.650762120863413e-06 = 1.0 field decay(t = 250.1): 0.0005431713534528096 / 0.0005431713534528096 = 1.0 field decay(t = 300.12): 0.015624179670709025 / 0.015624179670709025 = 1.0 field decay(t = 350.14): 0.10193257677430007 / 0.10193257677430007 = 1.0 field decay(t = 400.16): 0.23169084040266913 / 0.23169084040266913 = 1.0 field decay(t = 450.18): 0.23240467775624984 / 0.23240467775624984 = 1.0 field decay(t = 500.2): 0.10222488043703631 / 0.23240467775624984 = 0.4398572413600538 field decay(t = 550.22): 0.011956286799054389 / 0.23240467775624984 = 0.05144598170091204 field decay(t = 600.24): 0.0003757367393708341 / 0.23240467775624984 = 0.0016167348394119394 field decay(t = 650.26): 2.4029564352000626e-06 / 0.23240467775624984 = 1.0339535582499446e-05 field decay(t = 700.28): 4.4672075082859265e-09 / 0.23240467775624984 = 1.922167639401481e-08 run 0 finished at t = 700.28 (35014 timesteps) ----------- Initializing structure... time for choose_chunkdivision = 8.8819e-05 s Working in 1D dimensions. Computational cell is 0 x 0 x 100 with resolution 25 time for set_epsilon = 0.000489019 s ----------- field decay(t = 100.04): 5.34506997399672e-12 / 5.34506997399672e-12 = 1.0 field decay(t = 150.06): 1.016943844414871e-08 / 1.016943844414871e-08 = 1.0 field decay(t = 200.08): 4.649396296768624e-06 / 4.649396296768624e-06 = 1.0 field decay(t = 250.1): 0.0005265865745042597 / 0.0005265865745042597 = 1.0 field decay(t = 300.12): 0.012482616165581724 / 0.012482616165581724 = 1.0 field decay(t = 350.14): 0.07675630058831354 / 0.07675630058831354 = 1.0 field decay(t = 400.16): 0.18425111874447336 / 0.18425111874447336 = 1.0 field decay(t = 450.18): 0.20611275967726259 / 0.20611275967726259 = 1.0 field decay(t = 500.2): 0.17353357688171306 / 0.20611275967726259 = 0.8419351482821201 field decay(t = 550.22): 0.012903286335997464 / 0.20611275967726259 = 0.0626030448391541 field decay(t = 600.24): 0.00037570222684038757 / 0.20611275967726259 = 0.0018227994590372433 field decay(t = 650.26): 2.559656367429923e-06 / 0.20611275967726259 = 1.2418718624882361e-05 field decay(t = 700.28): 5.088653499594899e-09 / 0.20611275967726259 = 2.4688687432854045e-08 run 0 finished at t = 700.28 (35014 timesteps)
fig, ax = plt.subplots()
ax.semilogy(freqs, tflux[:, 0], "bo-", label=r"$\chi^{(3)}$=0.001")
ax.semilogy(freqs, tflux[:, 1], "ro-", label=r"$\chi^{(3)}$=0.01")
ax.semilogy(freqs, tflux[:, 2], "go-", label=r"$\chi^{(3)}$=0.1")
ax.semilogy(freqs, tflux[:, 3], "co-", label=r"$\chi^{(3)}$=1")
ax.set_xlabel("frequency")
ax.set_ylabel("transmitted power (a.u.)")
ax.set_xlim(0.2, 1.2)
ax.set_ylim(1e-15, 1e2)
ax.legend()
ax.grid(True)
plt.show()
/home/oskooi/.local/lib/python3.8/site-packages/numpy/core/getlimits.py:500: UserWarning: The value of the smallest subnormal for <class 'numpy.float64'> type is zero. setattr(self, word, getattr(machar, word).flat[0]) /home/oskooi/.local/lib/python3.8/site-packages/numpy/core/getlimits.py:89: UserWarning: The value of the smallest subnormal for <class 'numpy.float64'> type is zero. return self._float_to_str(self.smallest_subnormal)
For small values of $\chi^{(3)}$, we see a peak from our source at $\omega=\frac{1}{3}$ and another peak precisely at the third-harmonic frequency $3\omega=1$. As the $\chi^{(3)}$ gets larger, frequency-mixing within the peaks causes them to broaden, and finally for $\chi^{(3)}=1$ we start to see a noisy, broad-spectrum transmission due to the phenomenon of modulation instability. Notice also that at around $10^{-13}$ the data looks weird; this is probably due to our finite simulation time, imperfect absorbing boundaries, etcetera. We haven't attempted to analyze it in detail for this case.
In the second part of the tutorial, we investigate the dependence of the power at $\omega$ and $3\omega$ as a function of $\chi^{(3)}$ and the current amplitude. We could, of course, interpolate the flux spectrum above to get the desired frequencies, but it is easier just to add two flux regions to Meep and request exactly the desired frequency components.
logk = np.arange(-6.0, 0.2, 0.2)
first_order = np.zeros(len(logk))
third_order = np.zeros(len(logk))
for i, lk in enumerate(logk):
first_order[i], third_order[i] = third_harmonic_generation(
10**lk, flux_spectrum=False
)
----------- Initializing structure... time for choose_chunkdivision = 7.5685e-05 s Working in 1D dimensions. Computational cell is 0 x 0 x 100 with resolution 25 time for set_epsilon = 0.000831275 s ----------- field decay(t = 100.04): 5.345069974778501e-12 / 5.345069974778501e-12 = 1.0 field decay(t = 150.06): 1.0169446304565756e-08 / 1.0169446304565756e-08 = 1.0 field decay(t = 200.08): 4.6509134041875294e-06 / 4.6509134041875294e-06 = 1.0 field decay(t = 250.1): 0.0005451393864268183 / 0.0005451393864268183 = 1.0 field decay(t = 300.12): 0.01779607955736104 / 0.01779607955736104 = 1.0 field decay(t = 350.14): 0.13167020369297935 / 0.13167020369297935 = 1.0 field decay(t = 400.16): 0.24987915109109593 / 0.24987915109109593 = 1.0 field decay(t = 450.18): 0.2492659288871275 / 0.24987915109109593 = 0.9975459248949311 field decay(t = 500.2): 0.11135321997058951 / 0.24987915109109593 = 0.44562829465510145 field decay(t = 550.22): 0.012792973351210858 / 0.24987915109109593 = 0.051196641637969435 field decay(t = 600.24): 0.00037675013914112524 / 0.24987915109109593 = 0.0015077293863695624 field decay(t = 650.26): 2.4030142102446377e-06 / 0.24987915109109593 = 9.616705514453244e-06 field decay(t = 700.28): 4.46711423590339e-09 / 0.24987915109109593 = 1.7877098655080908e-08 run 0 finished at t = 700.28 (35014 timesteps) harmonics:, 1e-06, 1.0, 225.16443832253304, 2.515333788184986e-08 ----------- Initializing structure... time for choose_chunkdivision = 0.000133811 s Working in 1D dimensions. Computational cell is 0 x 0 x 100 with resolution 25 time for set_epsilon = 0.000910593 s ----------- field decay(t = 100.04): 5.345069974778501e-12 / 5.345069974778501e-12 = 1.0 field decay(t = 150.06): 1.0169446304561077e-08 / 1.0169446304561077e-08 = 1.0 field decay(t = 200.08): 4.650913403302946e-06 / 4.650913403302946e-06 = 1.0 field decay(t = 250.1): 0.0005451393753743672 / 0.0005451393753743672 = 1.0 field decay(t = 300.12): 0.017796064890280856 / 0.017796064890280856 = 1.0 field decay(t = 350.14): 0.13166946381281386 / 0.13166946381281386 = 1.0 field decay(t = 400.16): 0.24987671044152449 / 0.24987671044152449 = 1.0 field decay(t = 450.18): 0.24926351430698057 / 0.24987671044152449 = 0.9975460052541094 field decay(t = 500.2): 0.11135262363861735 / 0.24987671044152449 = 0.44563026078685236 field decay(t = 550.22): 0.012792966079088641 / 0.24987671044152449 = 0.0511971125939823 field decay(t = 600.24): 0.0003767501333900365 / 0.24987671044152449 = 0.0015077440899727334 field decay(t = 650.26): 2.403014206912613e-06 / 0.24987671044152449 = 9.616799431473868e-06 field decay(t = 700.28): 4.467114140778479e-09 / 0.24987671044152449 = 1.7877272887438072e-08 run 0 finished at t = 700.28 (35014 timesteps) harmonics:, 1.584893192461114e-06, 1.0, 225.1644324606108, 6.31859653887594e-08 ----------- Initializing structure... time for choose_chunkdivision = 8.0435e-05 s Working in 1D dimensions. Computational cell is 0 x 0 x 100 with resolution 25 time for set_epsilon = 0.000749297 s ----------- field decay(t = 100.04): 5.345069974778501e-12 / 5.345069974778501e-12 = 1.0 field decay(t = 150.06): 1.0169446304553853e-08 / 1.0169446304553853e-08 = 1.0 field decay(t = 200.08): 4.65091340190097e-06 / 4.65091340190097e-06 = 1.0 field decay(t = 250.1): 0.0005451393578573992 / 0.0005451393578573992 = 1.0 field decay(t = 300.12): 0.01779604164417494 / 0.01779604164417494 = 1.0 field decay(t = 350.14): 0.1316682910328865 / 0.1316682910328865 = 1.0 field decay(t = 400.16): 0.2498728412092159 / 0.2498728412092159 = 1.0 field decay(t = 450.18): 0.24925968639689838 / 0.2498728412092159 = 0.9975461326274986 field decay(t = 500.2): 0.11135167843143928 / 0.2498728412092159 = 0.44563337853194573 field decay(t = 550.22): 0.012792954553416579 / 0.2498728412092159 = 0.05119785924515571 field decay(t = 600.24): 0.0003767501242750534 / 0.2498728412092159 = 0.0015077674006179985 field decay(t = 650.26): 2.403014201626535e-06 / 0.2498728412092159 = 9.616948324586089e-06 field decay(t = 700.28): 4.4671139901279e-09 / 0.2498728412092159 = 1.787754911061996e-08 run 0 finished at t = 700.28 (35014 timesteps) harmonics:, 2.5118864315095823e-06, 1.0, 225.16442310949893, 1.5872094484029593e-07 ----------- Initializing structure... time for choose_chunkdivision = 8.1018e-05 s Working in 1D dimensions. Computational cell is 0 x 0 x 100 with resolution 25 time for set_epsilon = 0.000894799 s ----------- field decay(t = 100.04): 5.345069974778501e-12 / 5.345069974778501e-12 = 1.0 field decay(t = 150.06): 1.0169446304542361e-08 / 1.0169446304542361e-08 = 1.0 field decay(t = 200.08): 4.650913399679016e-06 / 4.650913399679016e-06 = 1.0 field decay(t = 250.1): 0.0005451393300948496 / 0.0005451393300948496 = 1.0 field decay(t = 300.12): 0.01779600480070021 / 0.01779600480070021 = 1.0 field decay(t = 350.14): 0.13166643192801916 / 0.13166643192801916 = 1.0 field decay(t = 400.16): 0.24986670621877566 / 0.24986670621877566 = 1.0 field decay(t = 450.18): 0.24925361691027295 / 0.24986670621877566 = 0.9975463345325971 field decay(t = 500.2): 0.11135018016629872 / 0.24986670621877566 = 0.44563832393421754 field decay(t = 550.22): 0.012792936286116477 / 0.24986670621877566 = 0.05119904320071908 field decay(t = 600.24): 0.0003767501098285919 / 0.24986670621877566 = 0.001507804363093981 field decay(t = 650.26): 2.403014193244504e-06 / 0.24986670621877566 = 9.617184416480433e-06 field decay(t = 700.28): 4.467113751670654e-09 / 0.24986670621877566 = 1.787798710469007e-08 run 0 finished at t = 700.28 (35014 timesteps) harmonics:, 3.981071705534978e-06, 1.0, 225.16440813678386, 3.9869373079802405e-07 ----------- Initializing structure... time for choose_chunkdivision = 7.4266e-05 s Working in 1D dimensions. Computational cell is 0 x 0 x 100 with resolution 25 time for set_epsilon = 0.000497459 s ----------- field decay(t = 100.04): 5.345069974778481e-12 / 5.345069974778481e-12 = 1.0 field decay(t = 150.06): 1.0169446304524124e-08 / 1.0169446304524124e-08 = 1.0 field decay(t = 200.08): 4.650913396157464e-06 / 4.650913396157464e-06 = 1.0 field decay(t = 250.1): 0.0005451392860940994 / 0.0005451392860940994 = 1.0 field decay(t = 300.12): 0.017795946405517928 / 0.017795946405517928 = 1.0 field decay(t = 350.14): 0.13166348450628076 / 0.13166348450628076 = 1.0 field decay(t = 400.16): 0.24985697620899291 / 0.24985697620899291 = 1.0 field decay(t = 450.18): 0.24924399074840453 / 0.24985697620899291 = 0.9975466546106135 field decay(t = 500.2): 0.11134780504186298 / 0.24985697620899291 = 0.4456461721874281 field decay(t = 550.22): 0.012792907333542469 / 0.24985697620899291 = 0.05120092113354417 field decay(t = 600.24): 0.00037675008693193487 / 0.24985697620899291 = 0.001507862988851679 field decay(t = 650.26): 2.403014179950846e-06 / 0.24985697620899291 = 9.617558878727662e-06 field decay(t = 700.28): 4.467113374285032e-09 / 0.24985697620899291 = 1.7878681804539708e-08 run 0 finished at t = 700.28 (35014 timesteps) harmonics:, 6.309573444801943e-06, 1.0, 225.16438402434432, 1.0014683697174824e-06 ----------- Initializing structure... time for choose_chunkdivision = 7.3466e-05 s Working in 1D dimensions. Computational cell is 0 x 0 x 100 with resolution 25 time for set_epsilon = 0.00130423 s ----------- field decay(t = 100.04): 5.345069974778434e-12 / 5.345069974778434e-12 = 1.0 field decay(t = 150.06): 1.0169446304495113e-08 / 1.0169446304495113e-08 = 1.0 field decay(t = 200.08): 4.6509133905761386e-06 / 4.6509133905761386e-06 = 1.0 field decay(t = 250.1): 0.0005451392163574408 / 0.0005451392163574408 = 1.0 field decay(t = 300.12): 0.017795853849839657 / 0.017795853849839657 = 1.0 field decay(t = 350.14): 0.13165881079960295 / 0.13165881079960295 = 1.0 field decay(t = 400.16): 0.24984153835108858 / 0.24984153835108858 = 1.0 field decay(t = 450.18): 0.24922871755702478 / 0.24984153835108858 = 0.9975471621007928 field decay(t = 500.2): 0.11134403938185904 / 0.24984153835108858 = 0.4456586367371521 field decay(t = 550.22): 0.012792861444657525 / 0.24984153835108858 = 0.051203901197087655 field decay(t = 600.24): 0.00037675005064174207 / 0.24984153835108858 = 0.0015079560153536836 field decay(t = 650.26): 2.403014158847725e-06 / 0.24984153835108858 = 9.618153068970066e-06 field decay(t = 700.28): 4.467112777638655e-09 / 0.24984153835108858 = 1.787978415086953e-08 run 0 finished at t = 700.28 (35014 timesteps) harmonics:, 1.0000000000000021e-05, 1.0, 225.16434484850555, 2.515516990517998e-06 ----------- Initializing structure... time for choose_chunkdivision = 7.2692e-05 s Working in 1D dimensions. Computational cell is 0 x 0 x 100 with resolution 25 time for set_epsilon = 0.00131563 s ----------- field decay(t = 100.04): 5.345069974778606e-12 / 5.345069974778606e-12 = 1.0 field decay(t = 150.06): 1.0169446304449148e-08 / 1.0169446304449148e-08 = 1.0 field decay(t = 200.08): 4.650913381730357e-06 / 4.650913381730357e-06 = 1.0 field decay(t = 250.1): 0.000545139105831829 / 0.000545139105831829 = 1.0 field decay(t = 300.12): 0.017795707145032064 / 0.017795707145032064 = 1.0 field decay(t = 350.14): 0.13165139755400873 / 0.13165139755400873 = 1.0 field decay(t = 400.16): 0.24981702876003392 / 0.24981702876003392 = 1.0 field decay(t = 450.18): 0.24920446914199645 / 0.24981702876003392 = 0.9975479669217191 field decay(t = 500.2): 0.11133806784509068 / 0.24981702876003392 = 0.4456784567397861 field decay(t = 550.22): 0.012792788710282837 / 0.24981702876003392 = 0.051208633669929574 field decay(t = 600.24): 0.00037674999312210095 / 0.24981702876003392 = 0.0015081037309269845 field decay(t = 650.26): 2.4030141253284044e-06 / 0.24981702876003392 = 9.619096573423189e-06 field decay(t = 700.28): 4.4671118357517325e-09 / 0.24981702876003392 = 1.7881534569217434e-08 run 0 finished at t = 700.28 (35014 timesteps) harmonics:, 1.5848931924611175e-05, 1.0, 225.16428034734233, 6.318402461852015e-06 ----------- Initializing structure... time for choose_chunkdivision = 7.5182e-05 s Working in 1D dimensions. Computational cell is 0 x 0 x 100 with resolution 25 time for set_epsilon = 0.000671585 s ----------- field decay(t = 100.04): 5.3450699747784145e-12 / 5.3450699747784145e-12 = 1.0 field decay(t = 150.06): 1.0169446304376287e-08 / 1.0169446304376287e-08 = 1.0 field decay(t = 200.08): 4.650913367710779e-06 / 4.650913367710779e-06 = 1.0 field decay(t = 250.1): 0.0005451389306594265 / 0.0005451389306594265 = 1.0 field decay(t = 300.12): 0.017795474598562566 / 0.017795474598562566 = 1.0 field decay(t = 350.14): 0.13163963349549604 / 0.13163963349549604 = 1.0 field decay(t = 400.16): 0.249778077761204 / 0.249778077761204 = 1.0 field decay(t = 450.18): 0.24916593257277922 / 0.249778077761204 = 0.9975492437370344 field decay(t = 500.2): 0.11132859514426717 / 0.249778077761204 = 0.44571003245008933 field decay(t = 550.22): 0.012792673420521304 / 0.249778077761204 = 0.05121615769960212 field decay(t = 600.24): 0.00037674990195060124 / 0.249778077761204 = 0.0015083385432679423 field decay(t = 650.26): 2.403014072018808e-06 / 0.249778077761204 = 9.620596385228682e-06 field decay(t = 700.28): 4.467110352557797e-09 / 0.249778077761204 = 1.7884317121010517e-08 run 0 finished at t = 700.28 (35014 timesteps) harmonics:, 2.5118864315095873e-05, 1.0, 225.1641720631706, 1.5869855336798398e-05 ----------- Initializing structure... time for choose_chunkdivision = 0.00018664 s Working in 1D dimensions. Computational cell is 0 x 0 x 100 with resolution 25 time for set_epsilon = 0.00048743 s ----------- field decay(t = 100.04): 5.345069974778471e-12 / 5.345069974778471e-12 = 1.0 field decay(t = 150.06): 1.0169446304260835e-08 / 1.0169446304260835e-08 = 1.0 field decay(t = 200.08): 4.650913345491174e-06 / 4.650913345491174e-06 = 1.0 field decay(t = 250.1): 0.0005451386530270641 / 0.0005451386530270641 = 1.0 field decay(t = 300.12): 0.017795105949300605 / 0.017795105949300605 = 1.0 field decay(t = 350.14): 0.131620951456545 / 0.131620951456545 = 1.0 field decay(t = 400.16): 0.2497160792387384 / 0.2497160792387384 = 1.0 field decay(t = 450.18): 0.24910459210960362 / 0.2497160792387384 = 0.997551270502889 field decay(t = 500.2): 0.11131356071984909 / 0.2497160792387384 = 0.44576048550493597 field decay(t = 550.22): 0.012792490664539549 / 0.2497160792387384 = 0.051228141589991186 field decay(t = 600.24): 0.0003767497574307566 / 0.2497160792387384 = 0.0015087124488710596 field decay(t = 650.26): 2.4030139870786908e-06 / 0.2497160792387384 = 9.622984608777694e-06 field decay(t = 700.28): 4.467108026452184e-09 / 0.2497160792387384 = 1.7888748053670394e-08 run 0 finished at t = 700.28 (35014 timesteps) harmonics:, 3.981071705534986e-05, 1.0, 225.16398523472878, 3.985811606969678e-05 ----------- Initializing structure... time for choose_chunkdivision = 7.3665e-05 s Working in 1D dimensions. Computational cell is 0 x 0 x 100 with resolution 25 time for set_epsilon = 0.000969299 s ----------- field decay(t = 100.04): 5.3450699747784185e-12 / 5.3450699747784185e-12 = 1.0 field decay(t = 150.06): 1.0169446304077684e-08 / 1.0169446304077684e-08 = 1.0 field decay(t = 200.08): 4.6509133102755144e-06 / 4.6509133102755144e-06 = 1.0 field decay(t = 250.1): 0.0005451382130023534 / 0.0005451382130023534 = 1.0 field decay(t = 300.12): 0.017794521458758514 / 0.017794521458758514 = 1.0 field decay(t = 350.14): 0.1315912490352465 / 0.1315912490352465 = 1.0 field decay(t = 400.16): 0.2496171544429268 / 0.2496171544429268 = 1.0 field decay(t = 450.18): 0.2490067133358019 / 0.2496171544429268 = 0.9975544905618077 field decay(t = 500.2): 0.11128967960573347 / 0.2496171544429268 = 0.44584147212998965 field decay(t = 550.22): 0.012792200930380799 / 0.2496171544429268 = 0.05124728289980425 field decay(t = 600.24): 0.00037674952832450226 / 0.2496171544429268 = 0.0015093094429559464 field decay(t = 650.26): 2.4030138513979544e-06 / 0.2496171544429268 = 9.626797712524147e-06 field decay(t = 700.28): 4.467104403481891e-09 / 0.2496171544429268 = 1.7895822959167907e-08 run 0 finished at t = 700.28 (35014 timesteps) harmonics:, 6.309573444801955e-05, 1.0, 225.16365094354785, 0.00010009819239837291 ----------- Initializing structure... time for choose_chunkdivision = 7.9811e-05 s Working in 1D dimensions. Computational cell is 0 x 0 x 100 with resolution 25 time for set_epsilon = 0.000878533 s ----------- field decay(t = 100.04): 5.345069974778442e-12 / 5.345069974778442e-12 = 1.0 field decay(t = 150.06): 1.0169446303787614e-08 / 1.0169446303787614e-08 = 1.0 field decay(t = 200.08): 4.650913254462442e-06 / 4.650913254462442e-06 = 1.0 field decay(t = 250.1): 0.0005451375155924395 / 0.0005451375155924395 = 1.0 field decay(t = 300.12): 0.017793594549247182 / 0.017793594549247182 = 1.0 field decay(t = 350.14): 0.131543940151223 / 0.131543940151223 = 1.0 field decay(t = 400.16): 0.24945871291377122 / 0.24945871291377122 = 1.0 field decay(t = 450.18): 0.2488499372366258 / 0.24945871291377122 = 0.9975596134926109 field decay(t = 500.2): 0.11125169749613016 / 0.24945871291377122 = 0.44597238635872305 field decay(t = 550.22): 0.01279174151811136 / 0.24945871291377122 = 0.05127799052876938 field decay(t = 600.24): 0.0003767491650682709 / 0.24945871291377122 = 0.0015102666115274128 field decay(t = 650.26): 2.403013633966078e-06 / 0.24945871291377122 = 9.632911217643907e-06 field decay(t = 700.28): 4.467098829144931e-09 / 0.24945871291377122 = 1.7907166989549264e-08 run 0 finished at t = 700.28 (35014 timesteps) harmonics:, 0.00010000000000000041, 1.0, 225.16302527038428, 0.00025135123845669517 ----------- Initializing structure... time for choose_chunkdivision = 7.78e-05 s Working in 1D dimensions. Computational cell is 0 x 0 x 100 with resolution 25 time for set_epsilon = 0.000898538 s ----------- field decay(t = 100.04): 5.3450699747783515e-12 / 5.3450699747783515e-12 = 1.0 field decay(t = 150.06): 1.0169446303327934e-08 / 1.0169446303327934e-08 = 1.0 field decay(t = 200.08): 4.650913166004646e-06 / 4.650913166004646e-06 = 1.0 field decay(t = 250.1): 0.0005451364102276239 / 0.0005451364102276239 = 1.0 field decay(t = 300.12): 0.017792124105148122 / 0.017792124105148122 = 1.0 field decay(t = 350.14): 0.13146837693935148 / 0.13146837693935148 = 1.0 field decay(t = 400.16): 0.24920348319619115 / 0.24920348319619115 = 1.0 field decay(t = 450.18): 0.24859736569326615 / 0.24920348319619115 = 0.9975677807743649 field decay(t = 500.2): 0.11119116730083108 / 0.24920348319619115 = 0.4461862485818197 field decay(t = 550.22): 0.012791012860002239 / 0.24920348319619115 = 0.05132758457445886 field decay(t = 600.24): 0.00037674858896835574 / 0.24920348319619115 = 0.0015118110876153035 field decay(t = 650.26): 2.403013284453902e-06 / 0.24920348319619115 = 9.642775669239241e-06 field decay(t = 700.28): 4.467090445070018e-09 / 0.24920348319619115 = 1.792547354385572e-08 run 0 finished at t = 700.28 (35014 timesteps) harmonics:, 0.00015848931924611207, 1.0, 225.16179312916677, 0.0006310268008721403 ----------- Initializing structure... time for choose_chunkdivision = 7.3666e-05 s Working in 1D dimensions. Computational cell is 0 x 0 x 100 with resolution 25 time for set_epsilon = 0.000957791 s ----------- field decay(t = 100.04): 5.345069974778269e-12 / 5.345069974778269e-12 = 1.0 field decay(t = 150.06): 1.0169446302599233e-08 / 1.0169446302599233e-08 = 1.0 field decay(t = 200.08): 4.650913025808378e-06 / 4.650913025808378e-06 = 1.0 field decay(t = 250.1): 0.0005451346582304672 / 0.0005451346582304672 = 1.0 field decay(t = 300.12): 0.017789790117491676 / 0.017789790117491676 = 1.0 field decay(t = 350.14): 0.1313471646264031 / 0.1313471646264031 = 1.0 field decay(t = 400.16): 0.24878880289605676 / 0.24878880289605676 = 1.0 field decay(t = 450.18): 0.24818694343126665 / 0.24878880289605676 = 0.9975808418313683 field decay(t = 500.2): 0.11109440471345641 / 0.24878880289605676 = 0.44654101559333975 field decay(t = 550.22): 0.012789856662544544 / 0.24878880289605676 = 0.05140848990654981 field decay(t = 600.24): 0.0003767476749381477 / 0.24878880289605676 = 0.0015143272950895294 field decay(t = 650.26): 2.403012722618603e-06 / 0.24878880289605676 = 9.658845955469203e-06 field decay(t = 700.28): 4.467078399956638e-09 / 0.24878880289605676 = 1.7955303244989568e-08 run 0 finished at t = 700.28 (35014 timesteps) harmonics:, 0.00025118864315095926, 1.0, 225.1592372361368, 0.0015836937554106964 ----------- Initializing structure... time for choose_chunkdivision = 0.000152275 s Working in 1D dimensions. Computational cell is 0 x 0 x 100 with resolution 25 time for set_epsilon = 0.000551973 s ----------- field decay(t = 100.04): 5.3450699747781285e-12 / 5.3450699747781285e-12 = 1.0 field decay(t = 150.06): 1.01694463014444e-08 / 1.01694463014444e-08 = 1.0 field decay(t = 200.08): 4.650912803612141e-06 / 4.650912803612141e-06 = 1.0 field decay(t = 250.1): 0.0005451318812208306 / 0.0005451318812208306 = 1.0 field decay(t = 300.12): 0.017786082245869233 / 0.017786082245869233 = 1.0 field decay(t = 350.14): 0.1311514603213696 / 0.1311514603213696 = 1.0 field decay(t = 400.16): 0.24810671668218798 / 0.24810671668218798 = 1.0 field decay(t = 450.18): 0.24751171235528827 / 0.24810671668218798 = 0.9976018209629453 field decay(t = 500.2): 0.11094283579553336 / 0.24810671668218798 = 0.4471577282514502 field decay(t = 550.22): 0.012788020821010097 / 0.24810671668218798 = 0.051542420906689514 field decay(t = 600.24): 0.00037674622379114167 / 0.24810671668218798 = 0.001518484581268851 field decay(t = 650.26): 2.403011830218934e-06 / 0.24810671668218798 = 9.685396116450443e-06 field decay(t = 700.28): 4.467062823838071e-09 / 0.24810671668218798 = 1.8004602549958978e-08 run 0 finished at t = 700.28 (35014 timesteps) harmonics:, 0.00039810717055349936, 1.0, 225.15367584280935, 0.003972440547314146 ----------- Initializing structure... time for choose_chunkdivision = 7.5129e-05 s Working in 1D dimensions. Computational cell is 0 x 0 x 100 with resolution 25 time for set_epsilon = 0.000496899 s ----------- field decay(t = 100.04): 5.345069974777999e-12 / 5.345069974777999e-12 = 1.0 field decay(t = 150.06): 1.0169446299614097e-08 / 1.0169446299614097e-08 = 1.0 field decay(t = 200.08): 4.650912451454404e-06 / 4.650912451454404e-06 = 1.0 field decay(t = 250.1): 0.0005451274792506687 / 0.0005451274792506687 = 1.0 field decay(t = 300.12): 0.017780183756957307 / 0.017780183756957307 = 1.0 field decay(t = 350.14): 0.13083246794663117 / 0.13083246794663117 = 1.0 field decay(t = 400.16): 0.24702535367001618 / 0.24702535367001618 = 1.0 field decay(t = 450.18): 0.24642377271214497 / 0.24702535367001618 = 0.9975646995381906 field decay(t = 500.2): 0.11075433688466904 / 0.24702535367001618 = 0.44835210329308134 field decay(t = 550.22): 0.01278510270744565 / 0.24702535367001618 = 0.05175623682953763 field decay(t = 600.24): 0.00037674391757787606 / 0.24702535367001618 = 0.0015251224701458854 field decay(t = 650.26): 2.4030104816223913e-06 / 0.24702535367001618 = 9.72778885212084e-06 field decay(t = 700.28): 4.467048183917948e-09 / 0.24702535367001618 = 1.808335912711682e-08 run 0 finished at t = 700.28 (35014 timesteps) harmonics:, 0.0006309573444801969, 1.0, 225.14108459473545, 0.009955033096505642 ----------- Initializing structure... time for choose_chunkdivision = 7.9957e-05 s Working in 1D dimensions. Computational cell is 0 x 0 x 100 with resolution 25 time for set_epsilon = 0.00116261 s ----------- field decay(t = 100.04): 5.345069974777698e-12 / 5.345069974777698e-12 = 1.0 field decay(t = 150.06): 1.0169446296713276e-08 / 1.0169446296713276e-08 = 1.0 field decay(t = 200.08): 4.650911893321058e-06 / 4.650911893321058e-06 = 1.0 field decay(t = 250.1): 0.0005451205008237796 / 0.0005451205008237796 = 1.0 field decay(t = 300.12): 0.017770780535768905 / 0.017770780535768905 = 1.0 field decay(t = 350.14): 0.13039636481724098 / 0.13039636481724098 = 1.0 field decay(t = 400.16): 0.24555052973215005 / 0.24555052973215005 = 1.0 field decay(t = 450.18): 0.24496560125506933 / 0.24555052973215005 = 0.9976178895736093 field decay(t = 500.2): 0.11044114777393672 / 0.24555052973215005 = 0.4497695358035166 field decay(t = 550.22): 0.012780456545355702 / 0.24555052973215005 = 0.052048173381245816 field decay(t = 600.24): 0.00037674024813005765 / 0.24555052973215005 = 0.0015342677066956898 field decay(t = 650.26): 2.4030087682352195e-06 / 0.24555052973215005 = 9.786208854269038e-06 field decay(t = 700.28): 4.467052624217102e-09 / 0.24555052973215005 = 1.819198935994895e-08 run 0 finished at t = 700.28 (35014 timesteps) harmonics:, 0.001000000000000006, 1.0, 225.1117121217118, 0.024907534309575777 ----------- Initializing structure... time for choose_chunkdivision = 8.1e-05 s Working in 1D dimensions. Computational cell is 0 x 0 x 100 with resolution 25 time for set_epsilon = 0.000806385 s ----------- field decay(t = 100.04): 5.3450699747772626e-12 / 5.3450699747772626e-12 = 1.0 field decay(t = 150.06): 1.0169446292115817e-08 / 1.0169446292115817e-08 = 1.0 field decay(t = 200.08): 4.65091100873642e-06 / 4.65091100873642e-06 = 1.0 field decay(t = 250.1): 0.0005451094363065589 / 0.0005451094363065589 = 1.0 field decay(t = 300.12): 0.01775574104989064 / 0.01775574104989064 = 1.0 field decay(t = 350.14): 0.1297186505527772 / 0.1297186505527772 = 1.0 field decay(t = 400.16): 0.24330568181730292 / 0.24330568181730292 = 1.0 field decay(t = 450.18): 0.24272872460013623 / 0.24330568181730292 = 0.9976286734742187 field decay(t = 500.2): 0.10991015097765587 / 0.24330568181730292 = 0.4517368857016125 field decay(t = 550.22): 0.012773039874014007 / 0.24330568181730292 = 0.05249791035955018 field decay(t = 600.24): 0.0003767344101864842 / 0.24330568181730292 = 0.0015483995580069202 field decay(t = 650.26): 2.403007840661162e-06 / 0.24330568181730292 = 9.876497016890748e-06 field decay(t = 700.28): 4.467122470311085e-09 / 0.24330568181730292 = 1.8360123927008928e-08 run 0 finished at t = 700.28 (35014 timesteps) harmonics:, 0.0015848931924611238, 1.0, 225.0417986570744, 0.06213865760610543 ----------- Initializing structure... time for choose_chunkdivision = 7.9428e-05 s Working in 1D dimensions. Computational cell is 0 x 0 x 100 with resolution 25 time for set_epsilon = 0.000509082 s ----------- field decay(t = 100.04): 5.3450699747765105e-12 / 5.3450699747765105e-12 = 1.0 field decay(t = 150.06): 1.0169446284829322e-08 / 1.0169446284829322e-08 = 1.0 field decay(t = 200.08): 4.65090960675786e-06 / 4.65090960675786e-06 = 1.0 field decay(t = 250.1): 0.0005450918890396788 / 0.0005450918890396788 = 1.0 field decay(t = 300.12): 0.0177315669710856 / 0.0177315669710856 = 1.0 field decay(t = 350.14): 0.12870448940754486 / 0.12870448940754486 = 1.0 field decay(t = 400.16): 0.24009442281853152 / 0.24009442281853152 = 1.0 field decay(t = 450.18): 0.23951896785692936 / 0.24009442281853152 = 0.9976032139570477 field decay(t = 500.2): 0.109172861261678 / 0.24009442281853152 = 0.45470802686738443 field decay(t = 550.22): 0.012761153864268465 / 0.24009442281853152 = 0.05315056349273726 field decay(t = 600.24): 0.00037672517750389907 / 0.24009442281853152 = 0.0015690709225204952 field decay(t = 650.26): 2.4030109102183914e-06 / 0.24009442281853152 = 1.0008607788589234e-05 field decay(t = 700.28): 4.467285139279417e-09 / 0.24009442281853152 = 1.8606367806618675e-08 run 0 finished at t = 700.28 (35014 timesteps) harmonics:, 0.0025118864315095977, 1.0, 224.87354381468265, 0.1541822732049199 ----------- Initializing structure... time for choose_chunkdivision = 7.4041e-05 s Working in 1D dimensions. Computational cell is 0 x 0 x 100 with resolution 25 time for set_epsilon = 0.00129713 s ----------- field decay(t = 100.04): 5.345069974775395e-12 / 5.345069974775395e-12 = 1.0 field decay(t = 150.06): 1.0169446273280971e-08 / 1.0169446273280971e-08 = 1.0 field decay(t = 200.08): 4.650907384754755e-06 / 4.650907384754755e-06 = 1.0 field decay(t = 250.1): 0.0005450640504066537 / 0.0005450640504066537 = 1.0 field decay(t = 300.12): 0.017700448238569356 / 0.017700448238569356 = 1.0 field decay(t = 350.14): 0.127238739165511 / 0.127238739165511 = 1.0 field decay(t = 400.16): 0.23607685977478424 / 0.23607685977478424 = 1.0 field decay(t = 450.18): 0.235521457134509 / 0.23607685977478424 = 0.9976473651809623 field decay(t = 500.2): 0.1080748102635521 / 0.23607685977478424 = 0.4577950179727685 field decay(t = 550.22): 0.012741991260680789 / 0.23607685977478424 = 0.05397391033088361 field decay(t = 600.24): 0.0003767107272766363 / 0.23607685977478424 = 0.0015957122084562453 field decay(t = 650.26): 2.403014548579255e-06 / 0.23607685977478424 = 1.0178949986337987e-05 field decay(t = 700.28): 4.467243285980297e-09 / 0.23607685977478424 = 1.8922834242381985e-08 run 0 finished at t = 700.28 (35014 timesteps) harmonics:, 0.003981071705535002, 1.0, 224.46789519742043, 0.37853579152133876 ----------- Initializing structure... time for choose_chunkdivision = 7.4749e-05 s Working in 1D dimensions. Computational cell is 0 x 0 x 100 with resolution 25 time for set_epsilon = 0.000481817 s ----------- field decay(t = 100.04): 5.345069974773556e-12 / 5.345069974773556e-12 = 1.0 field decay(t = 150.06): 1.0169446254978055e-08 / 1.0169446254978055e-08 = 1.0 field decay(t = 200.08): 4.650903863075193e-06 / 4.650903863075193e-06 = 1.0 field decay(t = 250.1): 0.0005450198586531935 / 0.0005450198586531935 = 1.0 field decay(t = 300.12): 0.017651930818716285 / 0.017651930818716285 = 1.0 field decay(t = 350.14): 0.12528210541303175 / 0.12528210541303175 = 1.0 field decay(t = 400.16): 0.23118021684180318 / 0.23118021684180318 = 1.0 field decay(t = 450.18): 0.2306509947916604 / 0.23118021684180318 = 0.9977107814095316 field decay(t = 500.2): 0.1065630103223593 / 0.23118021684180318 = 0.4609521168296181 field decay(t = 550.22): 0.012710818493179152 / 0.23118021684180318 = 0.054982293324333964 field decay(t = 600.24): 0.00037668761551044324 / 0.23118021684180318 = 0.0016294111176831836 field decay(t = 650.26): 2.4030055990318482e-06 / 0.23118021684180318 = 1.0394512263461658e-05 field decay(t = 700.28): 4.467125754656144e-09 / 0.23118021684180318 = 1.932313160564687e-08 run 0 finished at t = 700.28 (35014 timesteps) harmonics:, 0.0063095734448019814, 1.0, 223.49965201077922, 0.9098111290723069 ----------- Initializing structure... time for choose_chunkdivision = 7.397e-05 s Working in 1D dimensions. Computational cell is 0 x 0 x 100 with resolution 25 time for set_epsilon = 0.000488487 s ----------- field decay(t = 100.04): 5.345069974770661e-12 / 5.345069974770661e-12 = 1.0 field decay(t = 150.06): 1.0169446225969908e-08 / 1.0169446225969908e-08 = 1.0 field decay(t = 200.08): 4.650898281484014e-06 / 4.650898281484014e-06 = 1.0 field decay(t = 250.1): 0.0005449496426246324 / 0.0005449496426246324 = 1.0 field decay(t = 300.12): 0.017569355547691595 / 0.017569355547691595 = 1.0 field decay(t = 350.14): 0.12281472609090383 / 0.12281472609090383 = 1.0 field decay(t = 400.16): 0.22531439709826265 / 0.22531439709826265 = 1.0 field decay(t = 450.18): 0.22482927464433725 / 0.22531439709826265 = 0.9978469087631634 field decay(t = 500.2): 0.10798295720602 / 0.22531439709826265 = 0.4792545820271182 field decay(t = 550.22): 0.012669501060369235 / 0.22531439709826265 = 0.05623032182379316 field decay(t = 600.24): 0.00037665107710743793 / 0.22531439709826265 = 0.001671668929984866 field decay(t = 650.26): 2.4030053494996236e-06 / 0.22531439709826265 = 1.0665121183763684e-05 field decay(t = 700.28): 4.467048459398488e-09 / 0.22531439709826265 = 1.9825845649136873e-08 run 0 finished at t = 700.28 (35014 timesteps) harmonics:, 0.010000000000000082, 1.0, 221.25691828918846, 2.094806452639111 ----------- Initializing structure... time for choose_chunkdivision = 7.5854e-05 s Working in 1D dimensions. Computational cell is 0 x 0 x 100 with resolution 25 time for set_epsilon = 0.000920106 s ----------- field decay(t = 100.04): 5.345069974766093e-12 / 5.345069974766093e-12 = 1.0 field decay(t = 150.06): 1.0169446179995062e-08 / 1.0169446179995062e-08 = 1.0 field decay(t = 200.08): 4.6508894349934504e-06 / 4.6508894349934504e-06 = 1.0 field decay(t = 250.1): 0.0005448379145725761 / 0.0005448379145725761 = 1.0 field decay(t = 300.12): 0.0174518240541605 / 0.0174518240541605 = 1.0 field decay(t = 350.14): 0.1198104417750247 / 0.1198104417750247 = 1.0 field decay(t = 400.16): 0.2178144806969058 / 0.2178144806969058 = 1.0 field decay(t = 450.18): 0.2174095389093573 / 0.2178144806969058 = 0.9981408867479662 field decay(t = 500.2): 0.10689615163941552 / 0.2178144806969058 = 0.4907669650676905 field decay(t = 550.22): 0.012603441555172843 / 0.2178144806969058 = 0.057863194011930005 field decay(t = 600.24): 0.0003765927902040865 / 0.2178144806969058 = 0.0017289612196542828 field decay(t = 650.26): 2.4030010290977327e-06 / 0.2178144806969058 = 1.1032329078439774e-05 field decay(t = 700.28): 4.467142607822792e-09 / 0.2178144806969058 = 2.0508933077038762e-08 run 0 finished at t = 700.28 (35014 timesteps) harmonics:, 0.01584893192461127, 1.0, 216.41415446933203, 4.429179313929301 ----------- Initializing structure... time for choose_chunkdivision = 0.000131984 s Working in 1D dimensions. Computational cell is 0 x 0 x 100 with resolution 25 time for set_epsilon = 0.000912022 s ----------- field decay(t = 100.04): 5.345069974758849e-12 / 5.345069974758849e-12 = 1.0 field decay(t = 150.06): 1.016944610712986e-08 / 1.016944610712986e-08 = 1.0 field decay(t = 200.08): 4.650875413586146e-06 / 4.650875413586146e-06 = 1.0 field decay(t = 250.1): 0.0005446597284017753 / 0.0005446597284017753 = 1.0 field decay(t = 300.12): 0.01727707174154046 / 0.01727707174154046 = 1.0 field decay(t = 350.14): 0.11599910493080164 / 0.11599910493080164 = 1.0 field decay(t = 400.16): 0.20989110780110146 / 0.20989110780110146 = 1.0 field decay(t = 450.18): 0.2095734153022807 / 0.20989110780110146 = 0.9984863937203009 field decay(t = 500.2): 0.10413005947213903 / 0.20989110780110146 = 0.4961146785256645 field decay(t = 550.22): 0.012505896487971669 / 0.20989110780110146 = 0.059582783753862494 field decay(t = 600.24): 0.0003765002234598392 / 0.20989110780110146 = 0.001793788347701805 field decay(t = 650.26): 2.40299489953039e-06 / 0.20989110780110146 = 1.1448769434327506e-05 field decay(t = 700.28): 4.467105259129218e-09 / 0.20989110780110146 = 2.128296575271006e-08 run 0 finished at t = 700.28 (35014 timesteps) harmonics:, 0.025118864315096027, 1.0, 207.36551123280188, 7.9296875329420455 ----------- Initializing structure... time for choose_chunkdivision = 7.4101e-05 s Working in 1D dimensions. Computational cell is 0 x 0 x 100 with resolution 25 time for set_epsilon = 0.000909181 s ----------- field decay(t = 100.04): 5.345069974747352e-12 / 5.345069974747352e-12 = 1.0 field decay(t = 150.06): 1.0169445991646362e-08 / 1.0169445991646362e-08 = 1.0 field decay(t = 200.08): 4.6508531894836915e-06 / 4.6508531894836915e-06 = 1.0 field decay(t = 250.1): 0.0005443745521790347 / 0.0005443745521790347 = 1.0 field decay(t = 300.12): 0.0170313199261343 / 0.0170313199261343 = 1.0 field decay(t = 350.14): 0.111067032695486 / 0.111067032695486 = 1.0 field decay(t = 400.16): 0.20820740580453972 / 0.20820740580453972 = 1.0 field decay(t = 450.18): 0.20802140147634976 / 0.20820740580453972 = 0.9991066392308611 field decay(t = 500.2): 0.10095787089525556 / 0.20820740580453972 = 0.48489087362258604 field decay(t = 550.22): 0.012365563019300866 / 0.20820740580453972 = 0.05939060126857048 field decay(t = 600.24): 0.0003763526195198078 / 0.20820740580453972 = 0.0018075851724175405 field decay(t = 650.26): 2.4029889551311664e-06 / 0.20820740580453972 = 1.154132316209269e-05 field decay(t = 700.28): 4.467162646732303e-09 / 0.20820740580453972 = 2.145534943615777e-08 run 0 finished at t = 700.28 (35014 timesteps) harmonics:, 0.0398107170553501, 1.0, 193.6011512723667, 10.070607189474567 ----------- Initializing structure... time for choose_chunkdivision = 7.9907e-05 s Working in 1D dimensions. Computational cell is 0 x 0 x 100 with resolution 25 time for set_epsilon = 0.000907052 s ----------- field decay(t = 100.04): 5.345069974729155e-12 / 5.345069974729155e-12 = 1.0 field decay(t = 150.06): 1.0169445808617165e-08 / 1.0169445808617165e-08 = 1.0 field decay(t = 200.08): 4.650817962461843e-06 / 4.650817962461843e-06 = 1.0 field decay(t = 250.1): 0.0005439156829478632 / 0.0005439156829478632 = 1.0 field decay(t = 300.12): 0.016616296510880826 / 0.016616296510880826 = 1.0 field decay(t = 350.14): 0.10715059301879028 / 0.10715059301879028 = 1.0 field decay(t = 400.16): 0.20876708083835804 / 0.20876708083835804 = 1.0 field decay(t = 450.18): 0.20874885936981835 / 0.20876708083835804 = 0.9999127186696939 field decay(t = 500.2): 0.09911719000394087 / 0.20876708083835804 = 0.47477403815731023 field decay(t = 550.22): 0.01218124420974051 / 0.20876708083835804 = 0.05834849134654555 field decay(t = 600.24): 0.0003761164963470216 / 0.20876708083835804 = 0.0018016082556532803 field decay(t = 650.26): 2.4029780357721867e-06 / 0.20876708083835804 = 1.1510330202072131e-05 field decay(t = 700.28): 4.467009090915526e-09 / 0.20876708083835804 = 2.1397095140560952e-08 run 0 finished at t = 700.28 (35014 timesteps) harmonics:, 0.06309573444801994, 1.0, 173.54759042058004, 6.176862896998704 ----------- Initializing structure... time for choose_chunkdivision = 7.9734e-05 s Working in 1D dimensions. Computational cell is 0 x 0 x 100 with resolution 25 time for set_epsilon = 0.000541516 s ----------- field decay(t = 100.04): 5.345069974700306e-12 / 5.345069974700306e-12 = 1.0 field decay(t = 150.06): 1.0169445518535415e-08 / 1.0169445518535415e-08 = 1.0 field decay(t = 200.08): 4.650762120863385e-06 / 4.650762120863385e-06 = 1.0 field decay(t = 250.1): 0.0005431713534528077 / 0.0005431713534528077 = 1.0 field decay(t = 300.12): 0.015624179670708872 / 0.015624179670708872 = 1.0 field decay(t = 350.14): 0.1019325767743003 / 0.1019325767743003 = 1.0 field decay(t = 400.16): 0.23169084040267007 / 0.23169084040267007 = 1.0 field decay(t = 450.18): 0.23240467775625034 / 0.23240467775625034 = 1.0 field decay(t = 500.2): 0.10222488043703781 / 0.23240467775625034 = 0.4398572413600593 field decay(t = 550.22): 0.011956286799054538 / 0.23240467775625034 = 0.05144598170091257 field decay(t = 600.24): 0.00037573673937081094 / 0.23240467775625034 = 0.0016167348394118361 field decay(t = 650.26): 2.4029564352003662e-06 / 0.23240467775625034 = 1.033953558250073e-05 field decay(t = 700.28): 4.467207508276752e-09 / 0.23240467775625034 = 1.9221676393975294e-08 run 0 finished at t = 700.28 (35014 timesteps) harmonics:, 0.10000000000000102, 1.0, 143.29971061666146, 1.337907698289428 ----------- Initializing structure... time for choose_chunkdivision = 7.9257e-05 s Working in 1D dimensions. Computational cell is 0 x 0 x 100 with resolution 25 time for set_epsilon = 0.000487983 s ----------- field decay(t = 100.04): 5.345069974654548e-12 / 5.345069974654548e-12 = 1.0 field decay(t = 150.06): 1.0169445058786566e-08 / 1.0169445058786566e-08 = 1.0 field decay(t = 200.08): 4.650673591445848e-06 / 4.650673591445848e-06 = 1.0 field decay(t = 250.1): 0.0005419497773421297 / 0.0005419497773421297 = 1.0 field decay(t = 300.12): 0.01468960910671688 / 0.01468960910671688 = 1.0 field decay(t = 350.14): 0.09549018799286169 / 0.09549018799286169 = 1.0 field decay(t = 400.16): 0.2201120029446867 / 0.2201120029446867 = 1.0 field decay(t = 450.18): 0.22018394888915546 / 0.22018394888915546 = 1.0 field decay(t = 500.2): 0.09453094116088592 / 0.22018394888915546 = 0.4293271223347642 field decay(t = 550.22): 0.011701901150644354 / 0.22018394888915546 = 0.05314602272182565 field decay(t = 600.24): 0.00037512072334541086 / 0.22018394888915546 = 0.0017036697054345835 field decay(t = 650.26): 2.402942599168309e-06 / 0.22018394888915546 = 1.0913341373389545e-05 field decay(t = 700.28): 4.467083729548764e-09 / 0.22018394888915546 = 2.028796264253383e-08 run 0 finished at t = 700.28 (35014 timesteps) harmonics:, 0.15848931924611304, 1.0, 94.63060517366483, 0.3673651607719127 ----------- Initializing structure... time for choose_chunkdivision = 7.9286e-05 s Working in 1D dimensions. Computational cell is 0 x 0 x 100 with resolution 25 time for set_epsilon = 0.000483516 s ----------- field decay(t = 100.04): 5.345069974582112e-12 / 5.345069974582112e-12 = 1.0 field decay(t = 150.06): 1.0169444330133054e-08 / 1.0169444330133054e-08 = 1.0 field decay(t = 200.08): 4.650533215360899e-06 / 4.650533215360899e-06 = 1.0 field decay(t = 250.1): 0.0005399123569462003 / 0.0005399123569462003 = 1.0 field decay(t = 300.12): 0.014093771665812152 / 0.014093771665812152 = 1.0 field decay(t = 350.14): 0.11153132077038035 / 0.11153132077038035 = 1.0 field decay(t = 400.16): 0.20693495278636373 / 0.20693495278636373 = 1.0 field decay(t = 450.18): 0.2108832109073181 / 0.2108832109073181 = 1.0 field decay(t = 500.2): 0.1236622673624957 / 0.2108832109073181 = 0.5864016714770363 field decay(t = 550.22): 0.011423928445385981 / 0.2108832109073181 = 0.054171825230822805 field decay(t = 600.24): 0.0003741101509557853 / 0.2108832109073181 = 0.0017740158135215631 field decay(t = 650.26): 2.4027355648885305e-06 / 0.2108832109073181 = 1.1393678778651177e-05 field decay(t = 700.28): 4.466346571266359e-09 / 0.2108832109073181 = 2.1179242065075022e-08 run 0 finished at t = 700.28 (35014 timesteps) harmonics:, 0.2511886431509608, 1.0, 22.57026661243729, 0.413475056352828 ----------- Initializing structure... time for choose_chunkdivision = 7.8851e-05 s Working in 1D dimensions. Computational cell is 0 x 0 x 100 with resolution 25 time for set_epsilon = 0.000499921 s ----------- field decay(t = 100.04): 5.345069974467246e-12 / 5.345069974467246e-12 = 1.0 field decay(t = 150.06): 1.01694431752934e-08 / 1.01694431752934e-08 = 1.0 field decay(t = 200.08): 4.650310567520623e-06 / 4.650310567520623e-06 = 1.0 field decay(t = 250.1): 0.0005371041801518204 / 0.0005371041801518204 = 1.0 field decay(t = 300.12): 0.0133717723333235 / 0.0133717723333235 = 1.0 field decay(t = 350.14): 0.10093780191859889 / 0.10093780191859889 = 1.0 field decay(t = 400.16): 0.2156226336487837 / 0.2156226336487837 = 1.0 field decay(t = 450.18): 0.2158103663775883 / 0.2158103663775883 = 1.0 field decay(t = 500.2): 0.12022970754111681 / 0.2158103663775883 = 0.5571081202408938 field decay(t = 550.22): 0.011148924153637383 / 0.2158103663775883 = 0.05166074429497465 field decay(t = 600.24): 0.0003724922303640203 / 0.2158103663775883 = 0.001726016394005359 field decay(t = 650.26): 2.4037727963378503e-06 / 0.2158103663775883 = 1.1138356496425835e-05 field decay(t = 700.28): 4.458838350356822e-09 / 0.2158103663775883 = 2.066090904343077e-08 run 0 finished at t = 700.28 (35014 timesteps) harmonics:, 0.3981071705535018, 1.0, 9.726091154582134, 0.6671254223636301 ----------- Initializing structure... time for choose_chunkdivision = 7.8643e-05 s Working in 1D dimensions. Computational cell is 0 x 0 x 100 with resolution 25 time for set_epsilon = 0.000494216 s ----------- field decay(t = 100.04): 5.345069974285215e-12 / 5.345069974285215e-12 = 1.0 field decay(t = 150.06): 1.0169441344991794e-08 / 1.0169441344991794e-08 = 1.0 field decay(t = 200.08): 4.6499572759824876e-06 / 4.6499572759824876e-06 = 1.0 field decay(t = 250.1): 0.0005326910474981596 / 0.0005326910474981596 = 1.0 field decay(t = 300.12): 0.01292196509216716 / 0.01292196509216716 = 1.0 field decay(t = 350.14): 0.08319212231475127 / 0.08319212231475127 = 1.0 field decay(t = 400.16): 0.1911588547362798 / 0.1911588547362798 = 1.0 field decay(t = 450.18): 0.21391676622793007 / 0.21391676622793007 = 1.0 field decay(t = 500.2): 0.12912188463630342 / 0.21391676622793007 = 0.6036080617389429 field decay(t = 550.22): 0.011642738890991022 / 0.21391676622793007 = 0.05442649071548505 field decay(t = 600.24): 0.0003704269461133881 / 0.21391676622793007 = 0.0017316405471401672 field decay(t = 650.26): 2.407408397879035e-06 / 0.21391676622793007 = 1.1253949095854981e-05 field decay(t = 700.28): 4.514705562622321e-09 / 0.21391676622793007 = 2.110496359042687e-08 run 0 finished at t = 700.28 (35014 timesteps) harmonics:, 0.6309573444802007, 1.0, 22.671433283101095, 0.35168002847340085 ----------- Initializing structure... time for choose_chunkdivision = 7.6349e-05 s Working in 1D dimensions. Computational cell is 0 x 0 x 100 with resolution 25 time for set_epsilon = 0.000514033 s ----------- field decay(t = 100.04): 5.34506997399672e-12 / 5.34506997399672e-12 = 1.0 field decay(t = 150.06): 1.016943844414871e-08 / 1.016943844414871e-08 = 1.0 field decay(t = 200.08): 4.649396296768605e-06 / 4.649396296768605e-06 = 1.0 field decay(t = 250.1): 0.000526586574504257 / 0.000526586574504257 = 1.0 field decay(t = 300.12): 0.012482616165581575 / 0.012482616165581575 = 1.0 field decay(t = 350.14): 0.0767563005883142 / 0.0767563005883142 = 1.0 field decay(t = 400.16): 0.18425111874447933 / 0.18425111874447933 = 1.0 field decay(t = 450.18): 0.2061127596772555 / 0.2061127596772555 = 1.0 field decay(t = 500.2): 0.17353357688171156 / 0.2061127596772555 = 0.8419351482821418 field decay(t = 550.22): 0.012903286335997976 / 0.2061127596772555 = 0.06260304483915874 field decay(t = 600.24): 0.00037570222684037575 / 0.2061127596772555 = 0.0018227994590372485 field decay(t = 650.26): 2.5596563674296393e-06 / 0.2061127596772555 = 1.241871862488141e-05 field decay(t = 700.28): 5.0886534995752135e-09 / 0.2061127596772555 = 2.4688687432759386e-08 run 0 finished at t = 700.28 (35014 timesteps) harmonics:, 1.0000000000000122, 1.0, 9.438561367499654, 0.12758165549209102
We divide the transmitted power at all values of $\chi^{(3)}$ by the transmitted power at $\omega$ for the smallest $\chi^{(3)}$ of $10^{-6}$ which is 225.25726603587043. We then plot the fractional transmission at $\omega$ and $3\omega$ as a function of $\chi^{(3)}$ on a log-log scale.
input_flux = first_order[0]
fig, ax = plt.subplots()
ax.loglog(10**logk, first_order / input_flux, "ro-", label=r"$\omega$")
ax.loglog(10**logk, third_order / input_flux, "bo-", label=r"$3\omega$")
ax.loglog(10**logk, (10**logk) ** 2, "k-", label="quadratic line")
ax.set_xlabel(r"$\chi^{(3)}$")
ax.set_ylabel("transmission / incident power")
ax.legend()
ax.grid(True,'both')
plt.show()
As can be shown from coupled-mode theory or, equivalently, follows from Fermi's golden rule, the third-harmonic power must go as the square of $\chi^{(3)}$ as long as the nonlinearity is weak (i.e. in the first Born approximation limit, where the $\omega$ source is not depleted significantly). This is precisely what we see on the above graph, where the slope of the black line indicates an exact quadratic dependence, for comparison. Once the nonlinearity becomes strong enough, however, this approximation is no longer valid and the dependence is complicated.
Finally, we note that increasing the current amplitude by a factor of $F$ or the Kerr susceptibility $\chi^{(3)}$ by a factor $F^3$ should generate the same third-harmonic power in the weak nonlinearity approximation. And indeed, we see:
harmonics:, 0.001, 1.0, 225.2091048223644, 0.021498041565849526
harmonics:, 1e-06, 10.0, 22525.588597389557, 0.021791784143189268
which have third-harmonic powers differing by about 1% (last column).