Getting started¶

Let's start by importing all the pacakges. If error occurs, please make sure to run pip install -r requirements.txt, file which can be find at the root of this repository

In [1]:

import os 
import sys
import inspect

currentdir = os.getcwd()
parentdir = os.path.dirname(currentdir)
sys.path.insert(0, parentdir) 

import pyLLE
import numpy as np
from scipy import constants as cts
import pandas as pd
import pickle as pkl
import plotly.graph_objs as go
import plotly.graph_objs
import plotly.io as pio
import time
pio.renderers.default = "notebook"
pio.templates.default = "seaborn"

************************************************************
Julia solver:
/Users/greg/Documents/Lib/pyLLE_GitHub/pyLLE/ComputeLLE.jl
************************************************************

Here we will define the parameters. 2 dictionaries are needed, the resonator and the simulations. Novelty of the version 4 is that the pump power and the pump frequency can be lists to simulate a multi pump system:

In [2]:

res = dict(
        R=23e-6, 
        Qi=1e6, 
        Qc=1e6, 
        γ=3.2, 
        dispfile="./RW1000_H430.csv"
)

sim = dict(
    Pin=[140e-3], 
    f_pmp=[283e12],
    φ_pmp=[0], 
    δω=[None], 
    Tscan=0.7e6,
    μ_sim=[-220, 220],
    δω_init= 1e9 * 2 * np.pi,
    δω_end= -6.5e9 * 2 * np.pi,
    num_probe = 5000, 
)

The dispersion file which is available in the example folder is essentially a 2D column coma-separated-value (CSV) file that contain in the first column the azimuthal mode number and in the second the corresponding frequency of resonance. The value of the azimuthal mode number doesn't matter much given it is only to retrieve the effective refractive index and the LLE works in normalized mode number µ (normalized ot the pump). Therefore, as long as the mode number are monotonically increasing by 1, it is all good

Next step is to initialize the solver. Like in previous version, we do the same:

In [3]:

solver = pyLLE.LLEsolver(sim=sim, res=res,debug=False)

We need to analyze the linear dispersion of the resonator:

In [4]:

fig = solver.Analyze(plot=True)
fig.update_yaxes(range = [-50, 50], title = 'D<sub>int</sub> (GHz)')
fig.show()

Couple of New Stuff in v4¶

One first thing to notice compared to previous version is how all the parameters are stored in the solver object:

first, the solver has a custom print function that allows for easy display fo the parameters:

The simulation, resonator and dispersion parameters are stored as an attribute of the solver object, making it much easier to access all parameters:

In [5]:

solver.disp

Out[5]:

+------------------------------------------------------------------------------------------------------------------------+
|  Parameter      Description                           Values                                             Units         |
+------------------------------------------------------------------------------------------------------------------------+
|  ng             [1.760 ... 2.187]                     Group Index                                                      |
|  D              [-4108.943 ... -2129.830]             Dispersion                                         (ps/nm/km)    |
|  neff           [1.412 ... 1.936]                     Effctive Index                                                   |
|  freq           [148.354 ... 587.077] (THz)           Mode frequency                                     Hz            |
|  freq_sim       [67.530 ... 497.761] (THz)            Mode frequency to match the sim domain             Hz            |
|  Dint           [array for each pump]                 Integrated Dispersion at the pumps from file       Hz            |
|  Dint_sim       [array for each pump]                 Integrated Dispersion at the pumps from fit        Hz            |
|  Dint0          [-391064.774 ... -4964.929] (GHz)     Dint at the center of sim domain                   Hz            |
|  D1             6.172                                 angular repetition rate                            (x1e12  THz)  |
|  FSR            [982.263] (GHz)                       Free Spectra Range at the pumps                    Hz            |
|  FSR_center     982.263 (GHz)                         Free Spectra Range at the center of the domain     Hz            |
+------------------------------------------------------------------------------------------------------------------------+

In [6]:

solver.sim

Out[6]:

+---------------------------------------------------------------------------------------------------------+
|  Parameter          Description         Values                                      Units               |
+---------------------------------------------------------------------------------------------------------+
|  Pin                [140.000] (mW)      Pump power                                  W                   |
|  Tscan              700000.0            Simulation time length                      unit of round trip  |
|  f_pmp              [282.646] (THz)     Pump frequencies                            Hz                  |
|  φ_pmp              [0]                 Pump phase                                  rad                 |
|  f_center           282.646 (THz)       Center of the sim domain                    Hz                  |
|  δω_init            1.000 (GHz)         Start of the frequency sweep                x2π Hz              |
|  δω_end             -6.500 (GHz)        End of the frequency sweep                  x2π Hz              |
|  δω                 [None] (GHz)        Fixed detuning                              x2π Hz              |
|  μ_sim              [-220, 220]         Simulation mode domain                                          |
|  μ_fit              [None, None]        Fit mode domain                                                 |
|  μcenter            220                 Index of the center of the sim domain                           |
|  ind_pmp            [0]                 Pump index relative to center of domain                         |
|  ind_pump_sweep     0                   Pump index to sweep                                             |
|  num_probe          5000                Sampling number                                                 |
+---------------------------------------------------------------------------------------------------------+

In [7]:

solver.res

Out[7]:

+-----------------------------------------------------------------------------------+
|  Parameter     Description           Values                            Units      |
+-----------------------------------------------------------------------------------+
|  R             23.0                  Ring radius                       µm         |
|  Qi            1.000 x10^6           Intrinsic quality factor                     |
|  Qc            1.000 x10^6           Coupling quality factor                      |
|  γ             3.2                   Non-linear coefficient            W^-1 m^-1  |
|  dispfile      ./RW1000_H430.csv     mode/resonance frequency file                |
+-----------------------------------------------------------------------------------+

Each of teh parameters listed in the table can simply be accessed as attribute such as:

In [8]:

print(f"R = {solver.res.R}")
print(f"Pin = {solver.sim.Pin}")
print(f"D1 = {solver.disp.D1*1e-9/(2*np.pi)} x 2π GHz")

R = 2.3e-05
Pin = [0.14]
D1 = 982.2629817183881 x 2π GHz

Therfeore to plot stuff it is quite convenient as the integrated dispersion for instance can be retrieved through solver.disp.Dint

Simulation¶

We will run the simulation, now you can specify where the julia bin is present.

In [9]:

solver.Setup(verbose = False)
solver.SolveTemporal(bin = "/opt/bin/julia")
solver.RetrieveData()

----------------------------------------------------------------------
2023-09-16 11:49:03
Computing LLE [**************************************************] 100%

Simulation Time 00h:06min:18.6s
----------------------------------------------------------------------

Please note that you might have an error when you RetrieveData() because the hdf5 file took too long on the julia side to write and python is trying to open it. A quick workaround is to wait a bit (using time.sleep(XX) from the time package) or you can make a try/except loop. It will be fixed in further release hopefully

From here, all the solution components Are stored in the sol attribute of the solver object

In [10]:

solver.sol

Out[10]:

+-------------------------------------------------------------------------------------------------+
|  Parameter     Description                          Values                               Units  |
+-------------------------------------------------------------------------------------------------+
|  δfreq         [1.000 ... -6.499] (GHz)             Pump detuning                        Hz     |
|  time          [-0.095 ... 0.617] (μs)              Simualtion time                      s      |
|  Awg           [441 fast time x 5000 slow time]     E. field in time domain in wg,       V/m    |
|  Acav          [441 fast time x 5000 slow time]     E. field in time domain in cav.      V/m    |
|  Ewg           [441 spectra x 5000 slow time]       E. field in freq. domain in wg.      V/m    |
|  Ecav          [441 spectra x 5000 slow time]       E. field in freq. domain in cav.     V/m    |
|  Pcomb         [0.000 ... 0.000] (mW)               Intra-cavity comb power              W      |
|  Pwg           [139.242 ... 139.735] (mW)           In-waveguide power                   W      |
|  Pcav          [0.000 ... 0.265] (mW)               Intracavity power                    W      |
|  μ             [-220 ... 220]                       Frequency comb modes index                  |
|  freq          [66.548 ... 498.744] (THz)           Frequency comb frequencies           Hz     |
+-------------------------------------------------------------------------------------------------+

Retrieveing the intra cavity field power simply becomes:

In [11]:

fig = go.Figure()

tr = go.Scatter(y=solver.sol.Pcomb * 1e3)

fig.add_trace(tr)
fig.update_layout(
    xaxis_title="LLE sub-sampled step", yaxis_title="Intra-cavity Power (mW)"
)
fig.add_annotation(x=4500, y=4.9, ax=0, ay=-50, text="Single Soliton Step")
# fig.add_annotation(x=3151, y=7.6, ax=0, ay=-50, text="Two Solitons Step")
fig.add_annotation(x=1931, y=16, ax=0, ay=80, text="Modulation Instability")
fig.add_annotation(x=500, y=6, ax=-80, ay=-80, text="Primary Comb")
fig.show()

We can plot the spectra from here, either through the pre compile function or by processing it oursevelve from solver.sol.Ewg:

In [23]:

solver.PlotCombSpectra(3345).show()

Pretty spectra uh!?

Reusing a DKS as an initialize state¶

Sometime it is convenient to seed the simulation with a soliton, for convergence reason or inverse Fourier Transform. To this extend we have implemented the possibility to provide an input that is not zero but instead can be any random vector

In the following cell I am loading the solution from the previous simulation. I saved in a .csv file that is in this folder in which I also saved the detuning at which thte solution has been found. This is needed to be able to stay on the stable branch of the DKS

In [ ]:

df = pd.read_csv('./DKS_init.csv')
df["DKS_init"] = df.DKS_init.apply(lambda val: complex(val.strip('()')))
δω = df.det.values[0] * 2*np.pi
DKS_init = df.DKS_init.values

In [65]:

res = dict(
        R=23e-6, 
        Qi=1e6, 
        Qc=1e6, 
        γ=3.2, 
        dispfile="./RW1000_H430.csv"
)


sim = dict(
    Pin=[160e-3], # need to be the same length than fmp
    f_pmp=[283e12], # if only one pump f_pmp=[283e12]
    φ_pmp=[0], # need to be the same length than fmp
    δω=[None], # None defined which pump to sweep
    Tscan=0.2e6,
    μ_sim=[-220, 220],
    δω_end = δω, δω_init = δω, 
    DKS_init =  DKS_init, 
)
solver = pyLLE.LLEsolver(sim=sim, res=res,debug=False)

fig = solver.Analyze(plot=False)
solver.Setup(verbose = False)
solver.SolveTemporal()
solver.RetrieveData()

----------------------------------------------------------------------
2022-09-09 16:22:04
Computing LLE [************************************************* ] 99%

Simulation Time 00h:00min:25.1s
----------------------------------------------------------------------

Note that here I performed a much shorter simulation than in the detuning swept case. It is alright given the convergence of the solution is ok, we did not apply any perturbation to our system and the detuning is fixed. Essentially we are probing how stable the solution is

In the following cell, we will plot the temporal profile of the intracavity field Vs the LLE step and the resoantor angle:

In [66]:

fig = go.Figure()
x = np.arange(0, 5000)
y = np.linspace(-0.5, 0.5, solver.sol.Acav.shape[0])
xsteps = 10
ysteps = 2
tr = go.Heatmap(x= x[::xsteps], y = y[::ysteps],  z = np.abs(solver.sol.Acav[::ysteps, ::xsteps])**2)
fig.add_trace(tr)
fig.update_layout(xaxis_title = "LLE subsampled step", 
                 yaxis_title = "Resonator angle θ (x π)")
                 
fig.show()

We started with the soliton, and we kept the soliton all along. Great! Seems a bit underwhelming but think that you could potentially apply any perturbative approach through this solver feature (noise, etc...)

Freezing the soliton¶

In the previous simulation, we can notice that the soliton is actualy drifting with the LLE step, therefore in time. What it actualy means is the normalization of the LLE according to D1 is not the group rotation velocity coordinate of the soliton. It is well know that assymetry in the dispersion can lead to this kind of behavior. However, for many many reason (and I won't go in details here), it is super interesting to freeze the soliton envelope in time.

In the solver, this can be simply achieved by implementing a manual D1 value as followed:

In [63]:

df = pd.read_csv('./DKS_init.csv')
df["DKS_init"] = df.DKS_init.apply(lambda val: complex(val.strip('()')))
δω = df.det.values[0] * 2*np.pi
DKS_init = df.DKS_init.values

res = dict(
        R=23e-6, 
        Qi=1e6, 
        Qc=1e6, 
        γ=3.2, 
        dispfile="./RW1000_H430.csv"
)


sim = dict(
    Pin=[160e-3], # need to be the same length than fmp
    f_pmp=[283e12], # if only one pump f_pmp=[283e12]
    φ_pmp=[0], # need to be the same length than fmp
    δω=[None], # None defined which pump to sweep
    Tscan=0.2e6,
    μ_sim=[-220, 220],
    δω_end = δω, δω_init = δω, 
    DKS_init =  DKS_init, 
    D1_manual = 982.2557817183881 * 1e9 * 2*np.pi # linear is 982.2629817183881
)
solver = pyLLE.LLEsolver(sim=sim, res=res,debug=False)

fig = solver.Analyze(plot=False)
solver.Setup(verbose = False)
solver.SolveTemporal()
solver.RetrieveData()

----------------------------------------------------------------------
2022-09-09 16:18:08
Computing LLE [**************************************************] 100%

Simulation Time 00h:00min:25.9s
----------------------------------------------------------------------

In [64]:

fig = go.Figure()
x = np.arange(0, 5000)
y = np.linspace(-0.5, 0.5, solver.sol.Acav.shape[0])
xsteps = 10
ysteps = 2
tr = go.Heatmap(x= x[::xsteps], y = y[::ysteps],  z = np.abs(solver.sol.Acav[::ysteps, ::xsteps])**2)
fig.add_trace(tr)
fig.update_layout(xaxis_title = "LLE subsampled step", 
                 yaxis_title = "Resonator angle θ (x π)")
                 
fig.show()

Obviously, there is still a small drift, but to be honest I was lazy to actually compute correctly the shift of repetition rate and just did it manually. Long story short, you could figure it out and completely freezes the DKS in a coordinate system which moves at the same group velocity than the pulse

Multi Pump Simulations¶

In [96]:

df = pd.read_csv('./DKS_init.csv')
df["DKS_init"] = df.DKS_init.apply(lambda val: complex(val.strip('()')))
δω = df.det.values[0] * 2*np.pi
DKS_init = df.DKS_init.values

res = dict(
        R=23e-6, 
        Qi=1e6, 
        Qc=1e6, 
        γ=3.2, 
        dispfile="./RW1000_H430.csv"
)


sim = dict(
    Pin=[160e-3, 10e-3], # need to be the same length than fmp
    f_pmp=[283e12, 240e12], # if only one pump f_pmp=[283e12]
    φ_pmp=[0, 0], # need to be the same length than fmp
    δω=[None, 0], # None defined which pump to sweep
    Tscan=0.2e6,
    μ_sim=[-220, 220],
    δω_end = δω, δω_init = δω, 
    DKS_init =  DKS_init, 
    D1_manual = 982.2557817183881 * 1e9 * 2*np.pi # linear is 982.2629817183881
)
solver = pyLLE.LLEsolver(sim=sim, res=res,debug=False)

fig = solver.Analyze(plot=True)

In the case of multi pumps, the integratied dispersion will be displayed, relative to the group velocity (i.e. repetition rate) relative to the first pump. It is coherence with Moille et al, Nature Communications 2021:

In [97]:

fig = solver.Analyze(plot=True)
fig.update_yaxes(range = [-50, 50], title = 'D<sub>int</sub> (GHz)')
fig.show()

In [98]:

fig = solver.Analyze(plot=False)
solver.Setup(verbose = False)
solver.SolveTemporal()
solver.RetrieveData()

----------------------------------------------------------------------
2022-09-09 16:34:10
Computing LLE [************************************************* ] 98%

Simulation Time 00h:00min:24.9s
----------------------------------------------------------------------

In [102]:

fig = solver.PlotCombSpectra(4000)
fig.update_xaxes(range = [152, 460])
fig.show()

Couple of stuff to unpack here:

here I assumed that the second pump was small enough to consider it as a perturbation, which is semi-true. At high enough power, the repetition rate will be changed (see Miro Erkintalo papers, Hossein Taheri's among others)
We can clearly see new DW appearing. For complete explanation see our Nature comm paper (and soon a new one :) )