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, 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.
First, we'll load the necessary modules:
import meep as mp
import numpy as np
from matplotlib import pyplot as plt
%matplotlib notebook
Using MPI version 3.1, 1 processes
Next, we'll define some parameters of our simulation:
sz = 100 # size of cell in z direction
fcen = 1 / 3.0 # center frequency of source
df = fcen / 20.0 # frequency width of source
amp = 1 # amplitude of source
k = 10**-5 # Kerr susceptibility
dpml = 1.0 # PML thickness
Now, to define our cell, we'll do:
dimensions = 1
cell = mp.Vector3(0, 0, sz)
pml_layers = mp.PML(dpml)
resolution = 20
Note that this will only put PMLs at the $\pm z$ boundaries.
In this case, we're 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:
default_material = mp.Medium(index=1, chi3=k)
Now, our 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 our parameter amp
, above.
sources = mp.Source(mp.GaussianSource(fcen, fwidth=df), component=mp.Ex,
center=mp.Vector3(0, 0, -0.5*sz + dpml), amplitude=amp)
We'll want the frequency spectrum at the $+z$ end of the computational cell. 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
fmin = fcen / 2.0
fmax = fcen * 4
sim = mp.Simulation(cell_size=cell,
geometry=[],
sources=[sources],
boundary_layers=[pml_layers],
default_material=default_material,
resolution=resolution,
dimensions=dimensions)
trans = sim.add_flux(0.5 * (fmin + fmax), fmax - fmin, nfreq,
mp.FluxRegion(mp.Vector3(0, 0, 0.5*sz - dpml - 0.5)))
Finally, we'll run the sources, plus additional time for the field to decay at the flux plane, and output the flux spectrum:
sim.run(until_after_sources=mp.stop_when_fields_decayed(
50, mp.Ex, mp.Vector3(0, 0, 0.5*sz - dpml - 0.5), 1e-6))
----------- Initializing structure... field decay(t = 100.05000000000001): 4.1691008357827005e-12 / 4.1691008357827005e-12 = 1.0 field decay(t = 150.07500000000002): 1.0165420168812386e-08 / 1.0165420168812386e-08 = 1.0 field decay(t = 200.10000000000002): 4.650785262405232e-06 / 4.650785262405232e-06 = 1.0 field decay(t = 250.125): 0.0005454866366283557 / 0.0005454866366283557 = 1.0 field decay(t = 300.15000000000003): 0.017815669509667658 / 0.017815669509667658 = 1.0 field decay(t = 350.175): 0.13192124155368243 / 0.13192124155368243 = 1.0 field decay(t = 400.20000000000005): 0.2505510561658855 / 0.2505510561658855 = 1.0 field decay(t = 450.225): 0.2499503102652878 / 0.2505510561658855 = 0.9976023014638582 field decay(t = 500.25): 0.11167289094888443 / 0.2505510561658855 = 0.4457091207587955 field decay(t = 550.275): 0.012841517437730344 / 0.2505510561658855 = 0.05125309641172774 field decay(t = 600.3000000000001): 0.0003786349795018242 / 0.2505510561658855 = 0.0015112088741351644 field decay(t = 650.325): 2.4161065872412755e-06 / 0.2505510561658855 = 9.643170634417973e-06 field decay(t = 700.35): 4.490921803004586e-09 / 0.2505510561658855 = 1.7924178296144254e-08 run 0 finished at t = 700.35 (28014 timesteps)
In a linear calculation, we normalize the transmission against some reference spectrum, but in this case there is no obvious normalization so we will just plot the raw data. To do so, we'll pull the frequency points using get_flux_freqs()
and the corrensponding spectra using get_flux_freqs()
.
freqs = mp.get_flux_freqs(trans)
spectra = mp.get_fluxes(trans)
plt.figure(dpi=150)
plt.semilogy(freqs,spectra)
plt.grid(True)
plt.xlabel('Frequency')
plt.ylabel('Transmitted Power (a.u.)')
plt.show()
We next want to see what happens as we slowly increase our nonlinearity term ($\chi^{(3)}$). We'll wrap our routine in a function and parameterize it so that we can quickly loop over the various nonlinearities.
It is also interesting to have a more detailed look at 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 more flux regions to Meep and request exactly the desired frequency components. We'll add the additional fluxes to our function:
def run_chi3(k_pow,amp=1):
k = 10**k_pow
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)
sim = mp.Simulation(cell_size=cell,
geometry=[],
sources=[sources],
boundary_layers=[pml_layers],
default_material=default_material,
resolution=resolution,
dimensions=dimensions)
trans = sim.add_flux(0.5 * (fmin + fmax), fmax - fmin, nfreq,
mp.FluxRegion(mp.Vector3(0, 0, 0.5*sz - dpml - 0.5)))
# Single frequency point at omega
trans1 = sim.add_flux(fcen, 0, 1,
mp.FluxRegion(mp.Vector3(0, 0, 0.5*sz - dpml - 0.5)))
# Singel frequency point at 3omega
trans3 = sim.add_flux(3 * fcen, 0, 1,
mp.FluxRegion(mp.Vector3(0, 0, 0.5*sz - dpml - 0.5)))
sim.run(until_after_sources=mp.stop_when_fields_decayed(
50, mp.Ex, mp.Vector3(0, 0, 0.5*sz - dpml - 0.5), 1e-6))
omega_flux = mp.get_fluxes(trans1)
omega3_flux = mp.get_fluxes(trans3)
freqs = mp.get_flux_freqs(trans)
spectra = mp.get_fluxes(trans)
return freqs, spectra, omega_flux, omega3_flux
We'll now loop over various nonlinearities to see what effect this has on our frequency response.
k_pow = [-3,-2,-1,0]
freqs = []
spectra = []
for k_iter in k_pow:
freqs_iter, spectra_iter, omega_flux, omega3_flux = run_chi3(k_iter)
spectra.append(spectra_iter)
freqs = freqs_iter # Each iteration will simulate over the same frequencies, so just remember the last set.
----------- Initializing structure... field decay(t = 100.05000000000001): 4.1691008357817126e-12 / 4.1691008357817126e-12 = 1.0 field decay(t = 150.07500000000002): 1.0165420158581846e-08 / 1.0165420158581846e-08 = 1.0 field decay(t = 200.10000000000002): 4.650783221205379e-06 / 4.650783221205379e-06 = 1.0 field decay(t = 250.125): 0.0005454598673598967 / 0.0005454598673598967 = 1.0 field decay(t = 300.15000000000003): 0.017781404956060017 / 0.017781404956060017 = 1.0 field decay(t = 350.175): 0.13014055244440104 / 0.13014055244440104 = 1.0 field decay(t = 400.20000000000005): 0.2443206506136912 / 0.2443206506136912 = 1.0 field decay(t = 450.225): 0.2437445632391615 / 0.2443206506136912 = 0.9976420848050188 field decay(t = 500.25): 0.11044620451678894 / 0.2443206506136912 = 0.4520543156682302 field decay(t = 550.275): 0.012822890989720101 / 0.2443206506136912 = 0.05248386068681144 field decay(t = 600.3000000000001): 0.00037862010232358545 / 0.2443206506136912 = 0.0015496852246118257 field decay(t = 650.325): 2.41614906850782e-06 / 0.2443206506136912 = 9.8892543976077e-06 field decay(t = 700.35): 4.492194624953141e-09 / 0.2443206506136912 = 1.838647127727241e-08 run 0 finished at t = 700.35 (28014 timesteps) ----------- Initializing structure... field decay(t = 100.05000000000001): 4.169100835772216e-12 / 4.169100835772216e-12 = 1.0 field decay(t = 150.07500000000002): 1.0165420065576282e-08 / 1.0165420065576282e-08 = 1.0 field decay(t = 200.10000000000002): 4.650764664560235e-06 / 4.650764664560235e-06 = 1.0 field decay(t = 250.125): 0.0005452160206480376 / 0.0005452160206480376 = 1.0 field decay(t = 300.15000000000003): 0.017507545378757882 / 0.017507545378757882 = 1.0 field decay(t = 350.175): 0.11744036969268859 / 0.11744036969268859 = 1.0 field decay(t = 400.20000000000005): 0.20769602060422127 / 0.20769602060422127 = 1.0 field decay(t = 450.225): 0.2072148912173212 / 0.20769602060422127 = 0.9976834925122764 field decay(t = 500.25): 0.10101160503184105 / 0.20769602060422127 = 0.48634347802130234 field decay(t = 550.275): 0.012670991144119082 / 0.20769602060422127 = 0.06100738525108532 field decay(t = 600.3000000000001): 0.00037847776427704205 / 0.20769602060422127 = 0.001822267769868624 field decay(t = 650.325): 2.416057119618937e-06 / 0.20769602060422127 = 1.1632659656117805e-05 field decay(t = 700.35): 4.491338910111841e-09 / 0.20769602060422127 = 2.162457854053155e-08 run 0 finished at t = 700.35 (28014 timesteps) ----------- Initializing structure... field decay(t = 100.05000000000001): 4.169100835677209e-12 / 4.169100835677209e-12 = 1.0 field decay(t = 150.07500000000002): 1.0165419135520952e-08 / 1.0165419135520952e-08 = 1.0 field decay(t = 200.10000000000002): 4.650579070094068e-06 / 4.650579070094068e-06 = 1.0 field decay(t = 250.125): 0.0005427300142038431 / 0.0005427300142038431 = 1.0 field decay(t = 300.15000000000003): 0.015051531229532927 / 0.015051531229532927 = 1.0 field decay(t = 350.175): 0.10911319432715134 / 0.10911319432715134 = 1.0 field decay(t = 400.20000000000005): 0.22726682483704244 / 0.22726682483704244 = 1.0 field decay(t = 450.225): 0.22744269278874502 / 0.22744269278874502 = 1.0 field decay(t = 500.25): 0.10405323422921754 / 0.22744269278874502 = 0.45749209593585416 field decay(t = 550.275): 0.011388920462070745 / 0.22744269278874502 = 0.05007380242657031 field decay(t = 600.3000000000001): 0.00037704599306494387 / 0.22744269278874502 = 0.001657762614581575 field decay(t = 650.325): 2.4160426967236526e-06 / 0.22744269278874502 = 1.062264373983533e-05 field decay(t = 700.35): 4.491856927586553e-09 / 0.22744269278874502 = 1.9749400926055307e-08 run 0 finished at t = 700.35 (28014 timesteps) ----------- Initializing structure... field decay(t = 100.05000000000001): 4.16910083472712e-12 / 4.16910083472712e-12 = 1.0 field decay(t = 150.07500000000002): 1.0165409834939769e-08 / 1.0165409834939769e-08 = 1.0 field decay(t = 200.10000000000002): 4.648720329217878e-06 / 4.648720329217878e-06 = 1.0 field decay(t = 250.125): 0.0005182735319569587 / 0.0005182735319569587 = 1.0 field decay(t = 300.15000000000003): 0.013974697322011669 / 0.013974697322011669 = 1.0 field decay(t = 350.175): 0.07995289002032732 / 0.07995289002032732 = 1.0 field decay(t = 400.20000000000005): 0.17867229467477655 / 0.17867229467477655 = 1.0 field decay(t = 450.225): 0.23241654769547887 / 0.23241654769547887 = 1.0 field decay(t = 500.25): 0.18817882350661239 / 0.23241654769547887 = 0.8096618996043756 field decay(t = 550.275): 0.014319414613175686 / 0.23241654769547887 = 0.06161099437694745 field decay(t = 600.3000000000001): 0.0003754598154438489 / 0.23241654769547887 = 0.0016154607714756648 field decay(t = 650.325): 2.9312708100424914e-06 / 0.23241654769547887 = 1.261214332244172e-05 field decay(t = 700.35): 1.8599873539004834e-07 / 0.23241654769547887 = 8.002818096831515e-07 run 0 finished at t = 700.35 (28014 timesteps)
plt.figure(dpi=150)
plt.semilogy(freqs,np.array(spectra).T)
plt.grid(True)
plt.xlabel('Frequency')
plt.ylabel('Transmitted Power (a.u.)')
plt.legend(["$\chi^{{(3)}} = 10^{{{}}}$".format(i) for i in k_pow])
plt.show()
For small values of $\chi^{(3)}$, we see a peak from our source at $\omega=1/3$ and another peak precisely at the third-harmonic frequency $3\oemega=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.
Now, we can look specifically at our frequencies of interest. We'll run a quick simulation for a linear medium. We'll measure the PSD at $\omega$ and use this as a calibration factor for much larger nonlinearities.
_, _, omega_flux_cal, omega3_flux_cal = run_chi3(-16)
print("Omega: {}, 3Omega: {}".format(omega_flux_cal[0],omega3_flux_cal[0]))
----------- Initializing structure... field decay(t = 100.05000000000001): 4.1691008357827595e-12 / 4.1691008357827595e-12 = 1.0 field decay(t = 150.07500000000002): 1.016542016891584e-08 / 1.016542016891584e-08 = 1.0 field decay(t = 200.10000000000002): 4.650785283023354e-06 / 4.650785283023354e-06 = 1.0 field decay(t = 250.125): 0.0005454869069704206 / 0.0005454869069704206 = 1.0 field decay(t = 300.15000000000003): 0.017816014066636514 / 0.017816014066636514 = 1.0 field decay(t = 350.175): 0.13193941695549694 / 0.13193941695549694 = 1.0 field decay(t = 400.20000000000005): 0.2506140638300714 / 0.2506140638300714 = 1.0 field decay(t = 450.225): 0.2500128531691806 / 0.2506140638300714 = 0.9976010497906516 field decay(t = 500.25): 0.11168480654371428 / 0.2506140638300714 = 0.4456446092324733 field decay(t = 550.275): 0.012841704978003306 / 0.2506140638300714 = 0.05124095903376999 field decay(t = 600.3000000000001): 0.0003786351315011641 / 0.2506140638300714 = 0.001510829542901859 field decay(t = 650.325): 2.4161060812423408e-06 / 0.2506140638300714 = 9.640744195747047e-06 field decay(t = 700.35): 4.490904760943549e-09 / 0.2506140638300714 = 1.7919603921305078e-08 run 0 finished at t = 700.35 (28014 timesteps) Omega: 225.25726603587026, 3Omega: 5.026979907074879e-16
That is, the linear transmission is 225.25726603587026 at $\omega$, so we'll loop through several nonlinearities, divide by this value, and plot the fractional transmission at $\omega$ and $3\omega$ as a function of $\chi1^{(3)}$ on a log-log scale.
pts = np.linspace(-6,0,20)
_, _, omega_psd, omega3_psd = zip(*[run_chi3(k_iter) for k_iter in pts])
----------- Initializing structure... field decay(t = 100.05000000000001): 4.1691008357827595e-12 / 4.1691008357827595e-12 = 1.0 field decay(t = 150.07500000000002): 1.0165420168905492e-08 / 1.0165420168905492e-08 = 1.0 field decay(t = 200.10000000000002): 4.650785280961551e-06 / 4.650785280961551e-06 = 1.0 field decay(t = 250.125): 0.0005454868799362667 / 0.0005454868799362667 = 1.0 field decay(t = 300.15000000000003): 0.017815979612345514 / 0.017815979612345514 = 1.0 field decay(t = 350.175): 0.13193760003745314 / 0.13193760003745314 = 1.0 field decay(t = 400.20000000000005): 0.2506077676780644 / 0.2506077676780644 = 1.0 field decay(t = 450.225): 0.250006603482365 / 0.2506077676780644 = 0.9976011749305723 field decay(t = 500.25): 0.1116836154292776 / 0.2506077676780644 = 0.4456510524955018 field decay(t = 550.275): 0.012841686224530456 / 0.2506077676780644 = 0.051242171555620476 field decay(t = 600.3000000000001): 0.00037863511629851986 / 0.2506077676780644 = 0.0015108674396115361 field decay(t = 650.325): 2.416106131823971e-06 / 0.2506077676780644 = 9.640986607118051e-06 field decay(t = 700.35): 4.4909064662334595e-09 / 0.2506077676780644 = 1.7920060929645904e-08 run 0 finished at t = 700.35 (28014 timesteps) ----------- Initializing structure... field decay(t = 100.05000000000001): 4.1691008357827595e-12 / 4.1691008357827595e-12 = 1.0 field decay(t = 150.07500000000002): 1.0165420168894454e-08 / 1.0165420168894454e-08 = 1.0 field decay(t = 200.10000000000002): 4.650785278757183e-06 / 4.650785278757183e-06 = 1.0 field decay(t = 250.125): 0.0005454868510330065 / 0.0005454868510330065 = 1.0 field decay(t = 300.15000000000003): 0.0178159427756053 / 0.0178159427756053 = 1.0 field decay(t = 350.175): 0.13193565734819723 / 0.13193565734819723 = 1.0 field decay(t = 400.20000000000005): 0.2506010350871025 / 0.2506010350871025 = 1.0 field decay(t = 450.225): 0.24999992057177803 / 0.2506010350871025 = 0.9976013087291696 field decay(t = 500.25): 0.11168234185404839 / 0.2506010350871025 = 0.4456579431734209 field decay(t = 550.275): 0.012841666174342072 / 0.2506010350871025 = 0.05124346820785731 field decay(t = 600.3000000000001): 0.00037863510004548367 / 0.2506010350871025 = 0.0015109079653796314 field decay(t = 650.325): 2.416106185906883e-06 / 0.2506010350871025 = 9.641245835505452e-06 field decay(t = 700.35): 4.490908289112879e-09 / 0.2506010350871025 = 1.7920549639996314e-08 run 0 finished at t = 700.35 (28014 timesteps) ----------- Initializing structure... field decay(t = 100.05000000000001): 4.1691008357827595e-12 / 4.1691008357827595e-12 = 1.0 field decay(t = 150.07500000000002): 1.016542016887149e-08 / 1.016542016887149e-08 = 1.0 field decay(t = 200.10000000000002): 4.650785274196042e-06 / 4.650785274196042e-06 = 1.0 field decay(t = 250.125): 0.0005454867912281307 / 0.0005454867912281307 = 1.0 field decay(t = 300.15000000000003): 0.017815866554169706 / 0.017815866554169706 = 1.0 field decay(t = 350.175): 0.1319316371539862 / 0.1319316371539862 = 1.0 field decay(t = 400.20000000000005): 0.250587100703692 / 0.250587100703692 = 1.0 field decay(t = 450.225): 0.24998608899268374 / 0.250587100703692 = 0.997601585599097 field decay(t = 500.25): 0.11167970629203668 / 0.250587100703692 = 0.4456722073020547 field decay(t = 550.275): 0.012841624687286904 / 0.250587100703692 = 0.05124615214121316 field decay(t = 600.3000000000001): 0.00037863506641780544 / 0.250587100703692 = 0.0015109918481619068 field decay(t = 650.325): 2.416106297824614e-06 / 0.250587100703692 = 9.64178240236536e-06 field decay(t = 700.35): 4.49091206005567e-09 / 0.250587100703692 = 1.7921561195466212e-08 run 0 finished at t = 700.35 (28014 timesteps) ----------- Initializing structure... field decay(t = 100.05000000000001): 4.169100835782728e-12 / 4.169100835782728e-12 = 1.0 field decay(t = 150.07500000000002): 1.016542016882427e-08 / 1.016542016882427e-08 = 1.0 field decay(t = 200.10000000000002): 4.650785264758402e-06 / 4.650785264758402e-06 = 1.0 field decay(t = 250.125): 0.0005454866674834187 / 0.0005454866674834187 = 1.0 field decay(t = 300.15000000000003): 0.017815708836639688 / 0.017815708836639688 = 1.0 field decay(t = 350.175): 0.13192331666924295 / 0.13192331666924295 = 1.0 field decay(t = 400.20000000000005): 0.2505582526127686 / 0.2505582526127686 = 1.0 field decay(t = 450.225): 0.24995745365544636 / 0.2505582526127686 = 0.9976021585756716 field decay(t = 500.25): 0.11167425141378004 / 0.2505582526127686 = 0.4457017489915599 field decay(t = 550.275): 0.012841538842926424 / 0.2505582526127686 = 0.05125170976815797 field decay(t = 600.3000000000001): 0.000378634996846925 / 0.2505582526127686 = 0.0015111655389459305 field decay(t = 650.325): 2.4161065294750964e-06 / 0.2505582526127686 = 9.642893436079026e-06 field decay(t = 700.35): 4.49091985908172e-09 / 0.2505582526127686 = 1.7923655725769776e-08 run 0 finished at t = 700.35 (28014 timesteps) ----------- Initializing structure... field decay(t = 100.05000000000001): 4.1691008357828475e-12 / 4.1691008357828475e-12 = 1.0 field decay(t = 150.07500000000002): 1.0165420168726412e-08 / 1.0165420168726412e-08 = 1.0 field decay(t = 200.10000000000002): 4.650785245230697e-06 / 4.650785245230697e-06 = 1.0 field decay(t = 250.125): 0.0005454864114377906 / 0.0005454864114377906 = 1.0 field decay(t = 300.15000000000003): 0.017815382476511977 / 0.017815382476511977 = 1.0 field decay(t = 350.175): 0.13190609125009525 / 0.13190609125009525 = 1.0 field decay(t = 400.20000000000005): 0.2504984938450849 / 0.2504984938450849 = 1.0 field decay(t = 450.225): 0.24989813526348434 / 0.2504984938450849 = 0.9976033445455691 field decay(t = 500.25): 0.11166295794325522 / 0.2504984938450849 = 0.4457629913428168 field decay(t = 550.275): 0.012841361210895873 / 0.2504984938450849 = 0.05126322723056899 field decay(t = 600.3000000000001): 0.00037863485293490557 / 0.2504984938450849 = 0.0015115254671713265 field decay(t = 650.325): 2.4161070090188778e-06 / 0.2504984938450849 = 9.645195753205065e-06 field decay(t = 700.35): 4.490935979334019e-09 / 0.2504984938450849 = 1.792799593482321e-08 run 0 finished at t = 700.35 (28014 timesteps) ----------- Initializing structure... field decay(t = 100.05000000000001): 4.169100835782678e-12 / 4.169100835782678e-12 = 1.0 field decay(t = 150.07500000000002): 1.0165420168523849e-08 / 1.0165420168523849e-08 = 1.0 field decay(t = 200.10000000000002): 4.650785204825125e-06 / 4.650785204825125e-06 = 1.0 field decay(t = 250.125): 0.0005454858816409209 / 0.0005454858816409209 = 1.0 field decay(t = 300.15000000000003): 0.017814707103397807 / 0.017814707103397807 = 1.0 field decay(t = 350.175): 0.13187041021795567 / 0.13187041021795567 = 1.0 field decay(t = 400.20000000000005): 0.250374554159776 / 0.250374554159776 = 1.0 field decay(t = 450.225): 0.24977510744804346 / 0.250374554159776 = 0.9976058001830729 field decay(t = 500.25): 0.11163956208436386 / 0.250374554159776 = 0.4458902082083042 field decay(t = 550.275): 0.012840993630635223 / 0.250374554159776 = 0.051287135283087794 field decay(t = 600.3000000000001): 0.00037863455532979455 / 0.250374554159776 = 0.0015122725094826118 field decay(t = 650.325): 2.4161080020963407e-06 / 0.250374554159776 = 9.649974256387518e-06 field decay(t = 700.35): 4.490969255239367e-09 / 0.250374554159776 = 1.7937003503852332e-08 run 0 finished at t = 700.35 (28014 timesteps) ----------- Initializing structure... field decay(t = 100.05000000000001): 4.169100835782697e-12 / 4.169100835782697e-12 = 1.0 field decay(t = 150.07500000000002): 1.0165420168104837e-08 / 1.0165420168104837e-08 = 1.0 field decay(t = 200.10000000000002): 4.650785121220347e-06 / 4.650785121220347e-06 = 1.0 field decay(t = 250.125): 0.0005454847854047326 / 0.0005454847854047326 = 1.0 field decay(t = 300.15000000000003): 0.01781330928257498 / 0.01781330928257498 = 1.0 field decay(t = 350.175): 0.1317964137861031 / 0.1317964137861031 = 1.0 field decay(t = 400.20000000000005): 0.25011687030360397 / 0.25011687030360397 = 1.0 field decay(t = 450.225): 0.24951931317382292 / 0.25011687030360397 = 0.9976108883456933 field decay(t = 500.25): 0.11159103280876947 / 0.25011687030360397 = 0.4461555618911866 field decay(t = 550.275): 0.012840232906414287 / 0.25011687030360397 = 0.051336932574073034 field decay(t = 600.3000000000001): 0.0003786339402465762 / 0.25011687030360397 = 0.0015138280747994807 field decay(t = 650.325): 2.416110058754147e-06 / 0.25011687030360397 = 9.659924401826058e-06 field decay(t = 700.35): 4.491037714235971e-09 / 0.25011687030360397 = 1.795575687791284e-08 run 0 finished at t = 700.35 (28014 timesteps) ----------- Initializing structure... field decay(t = 100.05000000000001): 4.169100835782566e-12 / 4.169100835782566e-12 = 1.0 field decay(t = 150.07500000000002): 1.016542016723775e-08 / 1.016542016723775e-08 = 1.0 field decay(t = 200.10000000000002): 4.650784948230582e-06 / 4.650784948230582e-06 = 1.0 field decay(t = 250.125): 0.0005454825170836912 / 0.0005454825170836912 = 1.0 field decay(t = 300.15000000000003): 0.01781041537070784 / 0.01781041537070784 = 1.0 field decay(t = 350.175): 0.1316425937154995 / 0.1316425937154995 = 1.0 field decay(t = 400.20000000000005): 0.24957847104378408 / 0.24957847104378408 = 1.0 field decay(t = 450.225): 0.2489848350078133 / 0.24957847104378408 = 0.99762144533746 field decay(t = 500.25): 0.11149010814744109 / 0.24957847104378408 = 0.4467136435333084 field decay(t = 550.275): 0.012838658222915229 / 0.24957847104378408 = 0.05144136899798106 field decay(t = 600.3000000000001): 0.0003786326703923277 / 0.24957847104378408 = 0.0015170886687814647 field decay(t = 650.325): 2.4161143069600103e-06 / 0.24957847104378408 = 9.680780144438605e-06 field decay(t = 700.35): 4.4911771950517e-09 / 0.24957847104378408 = 1.7995050519657217e-08 run 0 finished at t = 700.35 (28014 timesteps) ----------- Initializing structure... field decay(t = 100.05000000000001): 4.1691008357823855e-12 / 4.1691008357823855e-12 = 1.0 field decay(t = 150.07500000000002): 1.0165420165443837e-08 / 1.0165420165443837e-08 = 1.0 field decay(t = 200.10000000000002): 4.650784590290659e-06 / 4.650784590290659e-06 = 1.0 field decay(t = 250.125): 0.0005454778233702294 / 0.0005454778233702294 = 1.0 field decay(t = 300.15000000000003): 0.017804420515894783 / 0.017804420515894783 = 1.0 field decay(t = 350.175): 0.13132132371111363 / 0.13132132371111363 = 1.0 field decay(t = 400.20000000000005): 0.24844276804284895 / 0.24844276804284895 = 1.0 field decay(t = 450.225): 0.24785728763909465 / 0.24844276804284895 = 0.9976433992892346 field decay(t = 500.25): 0.11127911965592859 / 0.24844276804284895 = 0.4479064556096729 field decay(t = 550.275): 0.01283539726004091 / 0.24844276804284895 = 0.05166339660902179 field decay(t = 600.3000000000001): 0.0003786300533999004 / 0.24844276804284895 = 0.0015240131817183668 field decay(t = 650.325): 2.416122931610749e-06 / 0.24844276804284895 = 9.725068476108912e-06 field decay(t = 700.35): 4.491452301556477e-09 / 0.24844276804284895 = 1.8078418369504868e-08 run 0 finished at t = 700.35 (28014 timesteps) ----------- Initializing structure... field decay(t = 100.05000000000001): 4.1691008357820155e-12 / 4.1691008357820155e-12 = 1.0 field decay(t = 150.07500000000002): 1.0165420161731733e-08 / 1.0165420161731733e-08 = 1.0 field decay(t = 200.10000000000002): 4.650783849662967e-06 / 4.650783849662967e-06 = 1.0 field decay(t = 250.125): 0.0005454681103843687 / 0.0005454681103843687 = 1.0 field decay(t = 300.15000000000003): 0.0177919867184528 / 0.0177919867184528 = 1.0 field decay(t = 350.175): 0.13064418737686334 / 0.13064418737686334 = 1.0 field decay(t = 400.20000000000005): 0.24608078239905923 / 0.24608078239905923 = 1.0 field decay(t = 450.225): 0.2454910539147474 / 0.24608078239905923 = 0.9976035167047076 field decay(t = 500.25): 0.11083352926803994 / 0.24608078239905923 = 0.45039489954280826 field decay(t = 550.275): 0.012828638329285625 / 0.24608078239905923 = 0.052131817057058695 field decay(t = 600.3000000000001): 0.00037862466733590196 / 0.24608078239905923 = 0.0015386194063781123 field decay(t = 650.325): 2.4161389249309716e-06 / 0.24608078239905923 = 9.818478718150437e-06 field decay(t = 700.35): 4.491929211746334e-09 / 0.24608078239905923 = 1.8253880566999964e-08 run 0 finished at t = 700.35 (28014 timesteps) ----------- Initializing structure... field decay(t = 100.05000000000001): 4.169100835781244e-12 / 4.169100835781244e-12 = 1.0 field decay(t = 150.07500000000002): 1.0165420154050911e-08 / 1.0165420154050911e-08 = 1.0 field decay(t = 200.10000000000002): 4.650782317199479e-06 / 4.650782317199479e-06 = 1.0 field decay(t = 250.125): 0.0005454480084043649 / 0.0005454480084043649 = 1.0 field decay(t = 300.15000000000003): 0.017766133980041062 / 0.017766133980041062 = 1.0 field decay(t = 350.175): 0.12943453172485406 / 0.12943453172485406 = 1.0 field decay(t = 400.20000000000005): 0.24162814253773343 / 0.24162814253773343 = 1.0 field decay(t = 450.225): 0.2410719556831018 / 0.24162814253773343 = 0.9976981702181286 field decay(t = 500.25): 0.10987489072588702 / 0.24162814253773343 = 0.45472720839514214 field decay(t = 550.275): 0.012814604477285866 / 0.24162814253773343 = 0.05303440378549736 field decay(t = 600.3000000000001): 0.00037861349073161515 / 0.24162814253773343 = 0.0015669262973889296 field decay(t = 650.325): 2.4161558408009105e-06 / 0.24162814253773343 = 9.999480256831406e-06 field decay(t = 700.35): 4.492861541918341e-09 / 0.24162814253773343 = 1.8594115299366343e-08 run 0 finished at t = 700.35 (28014 timesteps) ----------- Initializing structure... field decay(t = 100.05000000000001): 4.169100835779622e-12 / 4.169100835779622e-12 = 1.0 field decay(t = 150.07500000000002): 1.0165420138158315e-08 / 1.0165420138158315e-08 = 1.0 field decay(t = 200.10000000000002): 4.650779146309852e-06 / 4.650779146309852e-06 = 1.0 field decay(t = 250.125): 0.0005454063955002492 / 0.0005454063955002492 = 1.0 field decay(t = 300.15000000000003): 0.01771662173687749 / 0.01771662173687749 = 1.0 field decay(t = 350.175): 0.12689706943678725 / 0.12689706943678725 = 1.0 field decay(t = 400.20000000000005): 0.23320693569541642 / 0.23320693569541642 = 1.0 field decay(t = 450.225): 0.23267453504759023 / 0.23320693569541642 = 0.9977170462523398 field decay(t = 500.25): 0.10799807195922452 / 0.23320693569541642 = 0.46309974288362116 field decay(t = 550.275): 0.012785358594898285 / 0.23320693569541642 = 0.05482409241720325 field decay(t = 600.3000000000001): 0.000378589249816587 / 0.23320693569541642 = 0.0016234047614735157 field decay(t = 650.325): 2.416107138602139e-06 / 0.23320693569541642 = 1.0360357128304853e-05 field decay(t = 700.35): 4.493125804229867e-09 / 0.23320693569541642 = 1.9266690293028782e-08 run 0 finished at t = 700.35 (28014 timesteps) ----------- Initializing structure... field decay(t = 100.05000000000001): 4.169100835776254e-12 / 4.169100835776254e-12 = 1.0 field decay(t = 150.07500000000002): 1.0165420105274407e-08 / 1.0165420105274407e-08 = 1.0 field decay(t = 200.10000000000002): 4.6507725852541425e-06 / 4.6507725852541425e-06 = 1.0 field decay(t = 250.125): 0.0005453202108380395 / 0.0005453202108380395 = 1.0 field decay(t = 300.15000000000003): 0.017624098227291815 / 0.017624098227291815 = 1.0 field decay(t = 350.175): 0.12218039718656669 / 0.12218039718656669 = 1.0 field decay(t = 400.20000000000005): 0.21924792050420605 / 0.21924792050420605 = 1.0 field decay(t = 450.225): 0.21875381922070578 / 0.21924792050420605 = 0.9977463809811103 field decay(t = 500.25): 0.10436138938509092 / 0.21924792050420605 = 0.47599716861665214 field decay(t = 550.275): 0.012733297749038119 / 0.21924792050420605 = 0.05807716542877698 field decay(t = 600.3000000000001): 0.00037853854019597556 / 0.21924792050420605 = 0.0017265319521637774 field decay(t = 650.325): 2.4161208970651143e-06 / 0.21924792050420605 = 1.102004019700047e-05 field decay(t = 700.35): 4.491390315511021e-09 / 0.21924792050420605 = 2.0485440888935858e-08 run 0 finished at t = 700.35 (28014 timesteps) ----------- Initializing structure... field decay(t = 100.05000000000001): 4.169100835769312e-12 / 4.169100835769312e-12 = 1.0 field decay(t = 150.07500000000002): 1.0165420037232922e-08 / 1.0165420037232922e-08 = 1.0 field decay(t = 200.10000000000002): 4.650759009321587e-06 / 4.650759009321587e-06 = 1.0 field decay(t = 250.125): 0.000545141533396581 / 0.000545141533396581 = 1.0 field decay(t = 300.15000000000003): 0.017421291008439022 / 0.017421291008439022 = 1.0 field decay(t = 350.175): 0.11463570073897522 / 0.11463570073897522 = 1.0 field decay(t = 400.20000000000005): 0.20163815113463796 / 0.20163815113463796 = 1.0 field decay(t = 450.225): 0.20120253259282825 / 0.20163815113463796 = 0.9978396025783889 field decay(t = 500.25): 0.10109030457587416 / 0.20163815113463796 = 0.5013451274326259 field decay(t = 550.275): 0.01262521101057978 / 0.20163815113463796 = 0.06261320558404479 field decay(t = 600.3000000000001): 0.0003784354270167391 / 0.20163815113463796 = 0.0018768046864506805 field decay(t = 650.325): 2.4161436765026146e-06 / 0.20163815113463796 = 1.1982572062413455e-05 field decay(t = 700.35): 4.492170330890932e-09 / 0.20163815113463796 = 2.2278374928618628e-08 run 0 finished at t = 700.35 (28014 timesteps) ----------- Initializing structure... field decay(t = 100.05000000000001): 4.169100835754938e-12 / 4.169100835754938e-12 = 1.0 field decay(t = 150.07500000000002): 1.01654198964457e-08 / 1.01654198964457e-08 = 1.0 field decay(t = 200.10000000000002): 4.650730917976917e-06 / 4.650730917976917e-06 = 1.0 field decay(t = 250.125): 0.0005447703360772811 / 0.0005447703360772811 = 1.0 field decay(t = 300.15000000000003): 0.01702522296927303 / 0.01702522296927303 = 1.0 field decay(t = 350.175): 0.10518227880576582 / 0.10518227880576582 = 1.0 field decay(t = 400.20000000000005): 0.1772727477817735 / 0.1772727477817735 = 1.0 field decay(t = 450.225): 0.17689224598144007 / 0.1772727477817735 = 0.9978535798361865 field decay(t = 500.25): 0.0941896635654982 / 0.1772727477817735 = 0.5313262458223283 field decay(t = 550.275): 0.01239882951385902 / 0.1772727477817735 = 0.06994210711463807 field decay(t = 600.3000000000001): 0.00037822000681681363 / 0.1772727477817735 = 0.0021335485095679255 field decay(t = 650.325): 2.4160847007280536e-06 / 0.1772727477817735 = 1.3629194170907224e-05 field decay(t = 700.35): 4.493256393452313e-09 / 0.1772727477817735 = 2.5346571594769923e-08 run 0 finished at t = 700.35 (28014 timesteps) ----------- Initializing structure... field decay(t = 100.05000000000001): 4.169100835725167e-12 / 4.169100835725167e-12 = 1.0 field decay(t = 150.07500000000002): 1.016541960513757e-08 / 1.016541960513757e-08 = 1.0 field decay(t = 200.10000000000002): 4.650672789399802e-06 / 4.650672789399802e-06 = 1.0 field decay(t = 250.125): 0.000543995966140328 / 0.000543995966140328 = 1.0 field decay(t = 300.15000000000003): 0.016307143716002143 / 0.016307143716002143 = 1.0 field decay(t = 350.175): 0.09153211905255858 / 0.09153211905255858 = 1.0 field decay(t = 400.20000000000005): 0.22855260392168325 / 0.22855260392168325 = 1.0 field decay(t = 450.225): 0.22862376766772263 / 0.22862376766772263 = 1.0 field decay(t = 500.25): 0.0829714083179759 / 0.22862376766772263 = 0.36291680941311816 field decay(t = 550.275): 0.011966639931831521 / 0.22862376766772263 = 0.0523420642302755 field decay(t = 600.3000000000001): 0.0003777734295152949 / 0.22862376766772263 = 0.0016523803862087667 field decay(t = 650.325): 2.416110201923171e-06 / 0.22862376766772263 = 1.0568062221049122e-05 field decay(t = 700.35): 4.492316812744696e-09 / 0.22862376766772263 = 1.9649386669516103e-08 run 0 finished at t = 700.35 (28014 timesteps) ----------- Initializing structure... field decay(t = 100.05000000000001): 4.169100835663596e-12 / 4.169100835663596e-12 = 1.0 field decay(t = 150.07500000000002): 1.0165419002380455e-08 / 1.0165419002380455e-08 = 1.0 field decay(t = 200.10000000000002): 4.650552497483136e-06 / 4.650552497483136e-06 = 1.0 field decay(t = 250.125): 0.0005423672260522596 / 0.0005423672260522596 = 1.0 field decay(t = 300.15000000000003): 0.01462324952255367 / 0.01462324952255367 = 1.0 field decay(t = 350.175): 0.11286948557899079 / 0.11286948557899079 = 1.0 field decay(t = 400.20000000000005): 0.22340193165372965 / 0.22340193165372965 = 1.0 field decay(t = 450.225): 0.22423295751360775 / 0.22423295751360775 = 1.0 field decay(t = 500.25): 0.10875776760662008 / 0.22423295751360775 = 0.4850213314428591 field decay(t = 550.275): 0.011253173596244885 / 0.22423295751360775 = 0.050185190085458235 field decay(t = 600.3000000000001): 0.0003768391018113586 / 0.22423295751360775 = 0.0016805696450241478 field decay(t = 650.325): 2.416028682322177e-06 / 0.22423295751360775 = 1.0774636829091274e-05 field decay(t = 700.35): 4.491221029097323e-09 / 0.22423295751360775 = 2.002926366800817e-08 run 0 finished at t = 700.35 (28014 timesteps) ----------- Initializing structure... field decay(t = 100.05000000000001): 4.169100835536189e-12 / 4.169100835536189e-12 = 1.0 field decay(t = 150.07500000000002): 1.016541775519228e-08 / 1.016541775519228e-08 = 1.0 field decay(t = 200.10000000000002): 4.650303529004109e-06 / 4.650303529004109e-06 = 1.0 field decay(t = 250.125): 0.000538888737006022 / 0.000538888737006022 = 1.0 field decay(t = 300.15000000000003): 0.012940856099568963 / 0.012940856099568963 = 1.0 field decay(t = 350.175): 0.10610514124453076 / 0.10610514124453076 = 1.0 field decay(t = 400.20000000000005): 0.21497310015816543 / 0.21497310015816543 = 1.0 field decay(t = 450.225): 0.22343131384615714 / 0.22343131384615714 = 1.0 field decay(t = 500.25): 0.12060305444368463 / 0.22343131384615714 = 0.5397768664007654 field decay(t = 550.275): 0.010356579521407651 / 0.22343131384615714 = 0.046352408456670664 field decay(t = 600.3000000000001): 0.0003748549630386882 / 0.22343131384615714 = 0.0016777190116547107 field decay(t = 650.325): 2.4161257892828187e-06 / 0.22343131384615714 = 1.0813729497855586e-05 field decay(t = 700.35): 4.502116718465644e-09 / 0.22343131384615714 = 2.0149891440756423e-08 run 0 finished at t = 700.35 (28014 timesteps) ----------- Initializing structure... field decay(t = 100.05000000000001): 4.169100835272573e-12 / 4.169100835272573e-12 = 1.0 field decay(t = 150.07500000000002): 1.0165415174584912e-08 / 1.0165415174584912e-08 = 1.0 field decay(t = 200.10000000000002): 4.649788088402476e-06 / 4.649788088402476e-06 = 1.0 field decay(t = 250.125): 0.00053174759709517 / 0.00053174759709517 = 1.0 field decay(t = 300.15000000000003): 0.01247454839135486 / 0.01247454839135486 = 1.0 field decay(t = 350.175): 0.08976065034162566 / 0.08976065034162566 = 1.0 field decay(t = 400.20000000000005): 0.18870963593996914 / 0.18870963593996914 = 1.0 field decay(t = 450.225): 0.22826836798863853 / 0.22826836798863853 = 1.0 field decay(t = 500.25): 0.12497876557323505 / 0.22826836798863853 = 0.5475080348384299 field decay(t = 550.275): 0.009665581444509555 / 0.22826836798863853 = 0.042343061063066935 field decay(t = 600.3000000000001): 0.0003923542204648252 / 0.22826836798863853 = 0.0017188286923940053 field decay(t = 650.325): 2.7457665901575746e-06 / 0.22826836798863853 = 1.2028677535795227e-05 field decay(t = 700.35): 5.739133736098904e-09 / 0.22826836798863853 = 2.5142045683633898e-08 run 0 finished at t = 700.35 (28014 timesteps) ----------- Initializing structure... field decay(t = 100.05000000000001): 4.16910083472712e-12 / 4.16910083472712e-12 = 1.0 field decay(t = 150.07500000000002): 1.0165409834939769e-08 / 1.0165409834939769e-08 = 1.0 field decay(t = 200.10000000000002): 4.648720329217878e-06 / 4.648720329217878e-06 = 1.0 field decay(t = 250.125): 0.0005182735319569587 / 0.0005182735319569587 = 1.0 field decay(t = 300.15000000000003): 0.013974697322011669 / 0.013974697322011669 = 1.0 field decay(t = 350.175): 0.07995289002032732 / 0.07995289002032732 = 1.0 field decay(t = 400.20000000000005): 0.17867229467477655 / 0.17867229467477655 = 1.0 field decay(t = 450.225): 0.23241654769547887 / 0.23241654769547887 = 1.0 field decay(t = 500.25): 0.18817882350661239 / 0.23241654769547887 = 0.8096618996043756 field decay(t = 550.275): 0.014319414613175686 / 0.23241654769547887 = 0.06161099437694745 field decay(t = 600.3000000000001): 0.0003754598154438489 / 0.23241654769547887 = 0.0016154607714756648 field decay(t = 650.325): 2.9312708100424914e-06 / 0.23241654769547887 = 1.261214332244172e-05 field decay(t = 700.35): 1.8599873539004834e-07 / 0.23241654769547887 = 8.002818096831515e-07 run 0 finished at t = 700.35 (28014 timesteps)
quad = (10 ** pts) ** 2
plt.figure(dpi=150)
plt.loglog(10 ** pts,np.array(omega_psd)/omega_flux_cal[0],'o-',label='$\omega$')
plt.loglog(10 ** pts,np.array(omega3_psd)/omega_flux_cal[0],'o-',label='$3\omega$')
plt.loglog(10**pts,quad,'k',label='quadratic line')
plt.grid(True)
plt.xlabel('$\chi^{(3)}$')
plt.ylabel('Transmission/ Incident Power')
plt.legend()
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 gets 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:
_, _, omega_flux_1, omega3_flux_1 = run_chi3(-3,1)
_, _, omega_flux_2, omega3_flux_2 = run_chi3(-6,10)
print('-------------------------------')
print("Difference between powers: {}%".format(abs(omega3_flux_1[0]-omega3_flux_2[0])/omega3_flux_1[0] * 100))
----------- Initializing structure... field decay(t = 100.05000000000001): 4.1691008357817126e-12 / 4.1691008357817126e-12 = 1.0 field decay(t = 150.07500000000002): 1.0165420158581846e-08 / 1.0165420158581846e-08 = 1.0 field decay(t = 200.10000000000002): 4.650783221205379e-06 / 4.650783221205379e-06 = 1.0 field decay(t = 250.125): 0.0005454598673598967 / 0.0005454598673598967 = 1.0 field decay(t = 300.15000000000003): 0.017781404956060017 / 0.017781404956060017 = 1.0 field decay(t = 350.175): 0.13014055244440104 / 0.13014055244440104 = 1.0 field decay(t = 400.20000000000005): 0.2443206506136912 / 0.2443206506136912 = 1.0 field decay(t = 450.225): 0.2437445632391615 / 0.2443206506136912 = 0.9976420848050188 field decay(t = 500.25): 0.11044620451678894 / 0.2443206506136912 = 0.4520543156682302 field decay(t = 550.275): 0.012822890989720101 / 0.2443206506136912 = 0.05248386068681144 field decay(t = 600.3000000000001): 0.00037862010232358545 / 0.2443206506136912 = 0.0015496852246118257 field decay(t = 650.325): 2.41614906850782e-06 / 0.2443206506136912 = 9.8892543976077e-06 field decay(t = 700.35): 4.492194624953141e-09 / 0.2443206506136912 = 1.838647127727241e-08 run 0 finished at t = 700.35 (28014 timesteps) ----------- Initializing structure... field decay(t = 100.05000000000001): 4.169100835782652e-10 / 4.169100835782652e-10 = 1.0 field decay(t = 150.07500000000002): 1.0165420167882405e-06 / 1.0165420167882405e-06 = 1.0 field decay(t = 200.10000000000002): 0.0004650785076841832 / 0.0004650785076841832 = 1.0 field decay(t = 250.125): 0.05454842035005856 / 0.05454842035005856 = 1.0 field decay(t = 300.15000000000003): 1.7812567092193932 / 1.7812567092193932 = 1.0 field decay(t = 350.175): 13.175704416847648 / 13.175704416847648 = 1.0 field decay(t = 400.20000000000005): 24.99794167845129 / 24.99794167845129 = 1.0 field decay(t = 450.225): 24.9382864152439 / 24.99794167845129 = 0.9976135929919857 field decay(t = 500.25): 11.156520728608672 / 24.99794167845129 = 0.44629757410090326 field decay(t = 550.275): 1.2839829021813205 / 24.99794167845129 = 0.051363544994912066 field decay(t = 600.3000000000001): 0.03786336141250103 / 24.99794167845129 = 0.0015146591627237845 field decay(t = 650.325): 0.0002416111150361613 / 24.99794167845129 = 9.665240368347396e-06 field decay(t = 700.35): 4.491073802270905e-07 / 24.99794167845129 = 1.7965774382705666e-08 run 0 finished at t = 700.35 (28014 timesteps) ------------------------------- Difference between powers: 1.3663690082658835%
which have third-harmonic powers differing by about 1%.