pyLLE: a Fast and User Friendly Lugiato-Lefever Equation Solver

TL;DR

pyLLE is a user-friendly, cross-platform software tool developed in Python and Julia for efficiently simulating microresonator frequency combs using the Lugiato-Lefever equation.

Contribution

It introduces a new, accessible, and optimized solver for the Lugiato-Lefever equation, facilitating research in microresonator frequency combs.

Findings

Demonstrates ease of use and high performance compared to other languages.

Provides comprehensive documentation and examples.

Supports cross-platform operation.

Abstract

We present the development of pyLLE, a freely accessible and cross-platform Lugiato-Lefever equation solver programmed in Python and Julia and optimized for the simulation of microresonator frequency combs. Examples illustrating its operation, the simplicity of use, and performance against other programming language are presented. The documentation of the software can be found at https://gregmoille.github.io/pyLLE/

Tables

NIST Operating Unit	Physical Measurement Laboratory, Microsystems and Nanotechnology Division
Category	Lugiato Lefever Equation Solver.
Targeted Users	Scientist / Researcher / Engineer in Photonics/Metrology
Operating Systems	Windows/Linux/MacOS
Programming Language	Python$>$3.4, Julia 0.6.4
Inputs/Outputs	Input text file with resonator dispersion; output electric field data.
Documentation	https://usnistgov.github.io/pyLLE/
Accessibility	N/A.
Disclaimer	https://www.nist.gov/director/licensing

Table 2. Table 1: Performance comparison of the same resonator LLE simulation using three different software implementations.

Matlab 2018a	Pure Python	pyLLE
19 mins	45 min	11 mins

Equations

$$\frac{\partial E(z,\tau )}{\partial z}=-\frac{\alpha }{2}E+i\sum _{k>1}\frac{{\beta }_k}{k!}{\left(i\frac{\partial }{\partial \tau }\right)}^{k}E+i\gamma |E{|}^{2}E$$

$$E^{(m+1)}(0,\tau )=\theta \cdot E_{in}+(1-\theta )\cdot E^{(m)}\left(L,\tau \right)e^{i{\varphi }_0}$$

$$\frac{\partial E^{(m)}(z,\tau )}{\partial z}=-\frac{\alpha }{2}{E}^{(m)}+i\sum _{k>1}\frac{{\beta }_k}{k!}{\left(i\frac{\partial }{\partial \tau }\right)}^k{E}^{(m)}+i\gamma |E^{(m)}{|}^2{E}^{(m)}$$

$$t_R\frac{\partial E(t,\tau )}{\partial t}=\left(-\left(\frac{{\alpha }^{\prime }}{2}-i{\delta }_0\right)+i\sum _{k>1}\frac{{\beta }_k}{k!}{\left(i\frac{\partial }{\partial \tau }\right)}^k\right)E+\gamma |E{|}^{2}E+\sqrt{\theta }{E}_{in}$$

$$t_R\frac{\partial E(t,\tau )}{\partial t}=-\left(\frac{{\alpha }^{\prime }}{2}-i{\delta }_0\right)E+i\cdot {\mathrm{FT}}^{-1}\left[-t_R{D}_{int}(\omega )\cdot \mathrm{FT}\left[E(t,\tau )\right]\right]+\gamma |E{|}^{2}E+\sqrt{\theta }{E}_{in}$$

Full text

1 Introduction

The Lugiato-Lefever Equation (LLE), first developed to provide a description of spatial dissipative structures in optical systems [1], has recently made a significant impact in the integrated photonics community, where it has been adopted to help understand and predict Kerr-mediated nonlinear optical phenomena such as parametric frequency comb generation [2, 3] inside microresonators [4, 5, 6]. The LLE is essentially an application of the nonlinear Schrödinger equation (NLSE) [7] to a damped, driven Kerr nonlinear resonator, so that a periodic boundary condition is applied. Importantly, a slow-varying time envelope is stipulated, resulting in a mean-field solution in which the field does not vary within a round trip. This constraint, which differentiates the LLE from the more general Ikeda map [8], significantly simplifies calculations while still providing excellent physical representation for a wide variety of systems. In particular, simulations based on the LLE formalism have enabled modeling that quantitatively agrees with reported experimental results on microcomb generation [9, 10] (e.g., in terms of spectral bandwidth), and have also been central to theoretical studies that have provided better insight into novel nonlinear dynamics that can be supported by Kerr nonlinear microresonators [11, 12, 13].

The great potential of microresonator frequency combs (microcombs) in a wide variety of applications [3, 14, 15, 16, 17, 18, 19, 20, 21, 22] suggests the need for efficient and widely accessible computational tools to more rapidly further their development. Although LLE simulations are commonly performed by research groups working in the field, to our knowledge no free software package for solving this equation in an easy and fast way is currently available. Here, we introduce pyLLE, an open-source LLE solver for microcomb modeling. It combines the user-friendliness of the Python programming language and the computational power of the Julia programming language.

2 Software Specifications

3 Lugiato-Lefever Equation Model

The Lugiato-Lefever equation is derived from the nonlinear Schrödinger equation, which can be written, if considering only the third order nonlinearity and in the reference frame moving with the group velocity, as:

[TABLE]

with $E(z,\tau)$ being the complex electric field propagating along $z$ during the time $\tau$, $\alpha$ the loss per unit length, $\gamma=n_{2}\omega_{0}/(c_{0}A_{\rm eff})$ the effective nonlinear coefficient, $n_{2}$ is the Kerr nonlinear coefficient (i.e., the nonlinear refractive index), $\omega_{0}$ the angular frequency of nearest mode close to the pump frequency $\omega_{pmp}$, $A_{\rm eff}$ the effective mode area in the resonator at the pump frequency, $\beta_{\rm k}$ is the $k$th-order Taylor expansion coefficient of the dispersion [7], and $c_{0}$ the speed of light in vacuum.

Next, we assume a geometry such as the one described in Figure 1, where light injected into an optical waveguide is coupled into a traveling wave resonator (such as a ring or disk, for instance). Hence, a periodic condition for each round trip needs to be introduced, which leads to the coupled equations used to solve the Ikeda map [8]:

[TABLE]

where $E^{(m)}$ represents the field in the resonator for the $m$th round trip, $\phi_{0}$ the phase shift accumulated during a single round trip, $L=2\pi R$ the length of the resonator and $\theta$ the resonator-waveguide coupling coefficient.

We then make a key assumption, that of a slowly-varying field envelope, so that little non-linear physics happens during a single round trip. This assumption implies the possibility to average equations (2) and (3) over one round trip, leading to the LLE:

[TABLE]

with $t$ being the slow time and $\tau$ the fast time (i.e $-t_{R}/2\leq\tau\leq t_{R}/2$), $t_{R}$ the round trip time, $\alpha^{\prime}=\alpha L+\theta$ are the total round trip losses in the cavity and $\delta_{0}=\omega_{pmp}-\omega{{}_{0}}$ the detuning of the pump field with respect to the closest cold-cavity resonance $\omega_{0}$.

For systems with a complex dispersion profile, such as microring resonators, instead of using a Taylor expansion of the dispersion that results in a set of different $\beta$ coefficients (i.e., dispersion to different orders), it is often more convenient and accurate to work with the integrated dispersion $D_{int}=\omega_{\mu}-\left(\omega_{0}+D_{1}\mu\right)$, with $\mu$ an integer representing the mode number relative to the pump resonance (for which $\mu$=0), $\omega_{\mu}$ is the angular frequency of the $\mu$th mode, and $D_{1}$ is the angular free-spectral range evaluated at the pumped resonance (i.e., $D_{1}/2\pi$ is the free-spectral range in units of Hz). This implies that one has to work in the frequency domain. Considering the approximation of the slowly-varying envelope, one can consider the Fourier transform, $\mathrm{FT}$, and inverse Fourier transform, $\mathrm{FT}^{-1}$, relative to the fast time $\tau$. Hence, equation (4) becomes the one implemented in the package presented here:

[TABLE]

This approach is similar to the one found in Hansson and Wabnitz [8], for simulating the Kerr dynamics across a large number of modes using the Coupled Mode Theory (CMT) formalism. It is convenient indeed to work in the frequency domain, as it avoids the approximation of truncating the dispersion beyond a certain order, and the resonance frequencies can be found accurately using a Maxwell’s equation eigenvalue solver. Moreover, working in the frequency domain allows the use of well-developed computational libraries for executing the Fourier transform operation in one iteration of the split-step Fourier technique, and is the approach implemented in the pyLLE package.

4 pyLLE Description

4.1 Introduction to pyLLE

Several approaches for solving the complex field of a multi-mode resonator have been studied, including Coupled Mode Theory (CMT) and the LLE. However, the LLE provides a temporal approach which is well suited to study Kerr non-linear resonators, especially in the context of soliton microcombs. Our fundamental goal in the pyLLE package is to provide a fast, efficient, platform-independent and freely available tool for simulating the LLE, to enable researchers with interest in microcombs to simulate their behavior without requiring in-depth knowledge of the LLE formalism or its numerical implementation.

The bulk of the computation is done through the Julia back-end. The Julia language provides an efficient computational interface with BLAS (basic linear algebra subprograms), enabling reuse of the same Fast Fourier Transform plan (i.e plan_fft) throughout several iterations of a loop, resulting in an overall simulation that lasts only a few minutes. Python provides complementary functionality, with easy scripting, display of figures, and saving of results.

To assess the performance of the solver, we perform the same simulation of a Si3N4 ring resonator (described in more detail in the documentation of ref. [23]) using fixed hardware (a commercial laptop, the MacBook Pro 2017 3.1 GHz Intel Core i5, RAM 16 GB 2133 MHz LPDDR3) and three different software implementations, pyLLE, a pure Python version, and Matlab-based code [24]. We find that the simulation, in which the field behavior across 244 modes of the resonator is calculated, is an order of magnitude faster using pyLLE than the pure Python approach, and almost twice as fast as a Matlab solver (Table 1). We note that Matlab has a proprietary license, compared to the open-source pyLLE package.

As mentioned above, pyLLE is not only fast but highly user-friendly, thanks to the Python language which is becoming an increasingly valuable resource in the scientific toolbox. Below we describe a full example of how to use pyLLE, which consists ultimately of only a few lines of code. The complete example is accessible online [23].

4.2 Example of pyLLE use

In this section, we step through an example of how a simulation is set up and executed in pyLLE. We assume the pyLLE package has been installed on the user local machine following the instructions available in ref. [23].

4.2.1 Setting up pyLLE

Let’s start by importing the package:

⬇

import pyLLE

One has to define the resonator parameters.

⬇

res = {’R’: 23e-6, # ring radius in meters

       ’Qi’: 1e6,   # Intrinsic Q factor

       ’Qc’: 1e6,   # Coupling Q factor

       ’γ’:  1.55,  # Effective nonlinear coefficient at the pump frequency                            (units of W*-1*m*-1*)

       ’dispfile’: ’TestDispersion.txt’,

      }

The parameters in the code above are mostly self-explanatory and are chosen to fit the physical device investigated. The dispfile parameter should be a raw text file, without header, in a “comma separated value” (csv) format consisting of the integer azimuthal mode orders and corresponding resonance frequencies (in Hz) as the first and second column, respectively. Importantly, the parameter names are the ones commonly used in the LLE formalism, hence some follow the Greek alphabet. To improve the usability of pyLLE, one can define a parameter either by its Greek symbol or through its Latin name (e.g one can use either $\gamma$ or gamma).

The next step is to define the simulation parameters through the dictionary:

⬇

import numpy as np

sim = {’Pin’: 150e-3, # Input power in W

   ’Tscan’:   1e6,        # Simulation length in units of number of round trips

   ’f_pmp’:   191e12,     # Pump Frequency in Hz

   ’δω_init’: 2e9*2*np.pi,# Initial detuning of the pump in rad/s

   ’δω_end’: -8e9*2*np.pi,# End detuning of the pump in rad/s

   ’μ_sim’:  [-74,170],   # Azimuthal mode numbers to simulate in the LLE on the left and right side of the pump

   ’μ_fit’:  [-71,180],   # Azimuthal mode number range to fit the dispersion on the left and right side of the pump

   }

In the parameters described above, one specifies a simulation with pump power Pin in the input coupling bus waveguide of 150 mW, with a linear detuning ramp of the pump from $\delta\omega$_init to $\delta\omega$_end relative to the pump mode angular frequency, i.e. the mode closest to the defined pump frequency f_pmp. The simulation length Tscan is in units of number of round trips, as is most convenient in the LLE formalism. Two parameter ranges for the simulation frequency range have to be defined, μ_fit, which determines the fit window for the raw data found in dispfile, and μ_sim which is the number of modes simulated in the LLE. μ_sim can differ from μ_fit (i.e., the simulation can be performed over a truncated range or can be expanded over a larger range through extrapolation).

Once the resonator and simulation parameters are set, one can easily instantiate the pyLLE class:

⬇

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

4.2.2 Dispersion Analysis

To plot and retrieve all the dispersion data, the method solver.Analyze has to be called, resulting in the plot of the integrated dispersion $D_{int}$, and returning the handle of the figure (and axes if using matplotlib).

⬇

_ = solver.Analyze(plot=True, plottype=’all’)

Here, one extrapolates (orange curve - fit dispersion) outside of the spectral window containing the raw dispersion data (green circles), to enable a simulation over a broader spectral range (blue dashed line). One has to be careful about using this extrapolation feature, to make sure that the extrapolated dispersion is a reasonable physical representation of the dispersion. For example, extrapolation can result in ripples in the integrated dispersion, creating zero-crossings that are artifacts that will influence the LLE simulation results.

4.2.3 Temporal Simulation

The core function of pyLLE is to solve the full temporal Lugiato-Lefever Equation as described in eq. (5). The solver implemented in this package is based on a Julia core called from python. To interface both languages easily, a .hdf5 file is created in a temporary location with all the necessary data to solve the LLE.

We first set up the solver:

⬇

solver.Setup()

After which we call the temporal solver :

⬇

solver.SolveTemporal()

A progress bar has been implemented which displays the state of the current simulation, making it convenient to assess the time a simulation will take and when it is complete. To process the data through python, one has to retrieve the Julia data using the method:

⬇

solver.RetrieveData()

The pyLLE class includes methods to ease the plotting process. The solver.PlotCombPower method (fig.3) provides an overview of the microcomb behavior, plotting the spectra, time envelope, and comb power inside the resonator for the different steps of the LLE, which are sub-sampled for faster display.

⬇ solver.PlotCombPower()

To obtain a better idea of the microcomb behavior for a given laser-cavity mode detuning, we can plot the spectrum and the temporal profile of the electric field in the resonator at a fixed LLE step:

⬇ ind = 570

_ = solver.PlotCombSpectra(ind)

The temporal profile can also be retrieved, which can be interesting, for example, in the case when we are on a soliton state for the given detuning, and by comparing the position of the soliton relative to the origin of time which is linked to the soliton drift:

⬇ ind = 570

_ = solver.PlotSolitonTime(ind)

4.2.4 Steady-state Simulation

Instead of solving the full temporal LLE in eq (5), we can instead consider solutions to the steady-state LLE, i.e., considering the right-hand side of eq (5) equal to zero. Steady-state solutions to the LLE indicate what types of states are supported by the system, but do not necessarily give an indication of how to access such states in practice.

Newton’s method can be used to find steady-state solutions and is implemented in the method solver.SolveSteadyState. Although it gives fast results, the accuracy of such a solver remains questionable compared to a full temporal resolution of the equation (in particular, the convergence of the Newton’s method is not always assured).

To run the steady-state solver, one has to set it up in the same way as in the temporal method (see sections 4.2.1 to 4.2.2 ), i.e. defining the two dictionaries res and sim as in the example above, and calculating the dispersion profile. Here it is not necessary to use the method solver.Setup, as the Newton’s method places relatively low computational demands and thus everything is solved with Python.

However, the main difference with respect to the temporal solver is the introduction of a fixed detuning $\delta\omega$ of the pump relative to the closest mode, so that:

⬇

solver.sim[’δω’] = -5e92np.pi # this value corresponds to the detuning at the end of the soliton step in the temporal simulation

The results of the steady-state Newton’s method solver, with the function call indicated below, produces the spectrum shown in Fig. 6 and returns the figure handle:

⬇

steady_fig = solver.SolveSteadyState()

5 pyLLE Module Description

5.1 _llesolver.py

5.1.1 MyLogger(self, fname)

Helper class for logging the results

5.1.2 Latexify(self, **kwarg)

Class that handles saving the figures in a nice way compatible with the column/page size of different latex templates.

• •

Input ([] = optional):

• –

figname = name to save the figure (without extension)

• –

fig = matplotlib handle to the figure

• –

[fig_width]: default = ‘1column’

• –

[frmt]: default = ‘pdf’

• –

[fig_height] : default = 6.5

• –

[font_size] : default = 8

5.1.3 LLEsolver(self, **kwargs)

Class to solve the LLE Initialization input ([] = optional):

• •

res <dict>

• –

Qi <float>: intrinsic Q of the resonator

• –

Qc <float>: coupling Q of the resonator

• –

R <float>: ring radius

• –

gamma <float>: Effective non-linearity of the waveguide comprising the resonator.

• –

dispfile <str> : str pointing to a .csv file where the azimuthal mode orders and corresponding resonance frequencies are saved

• •

sim <dict>

• –

Tscan <float>: length of the simulation (in units of round trip)

• –

mu_fit <list>: number of modes to fit

• –

mu_sim <list>: number of mode that were simulated

• –

domega_init <float>: initial detuning of the pump

• –

domega_end <float>: final detuning of the pump

• –

[domga_stop] <float>: where to stop the scan in detuning but keep doing the simulation

• •

debug <bool>: Save a trace in a log file in the working directory of the different actions pyLLE perform (default = True)

LLEsolver.Analyze(self, plot=False, f=None, ax=None, label=None, plottype=‘all’, zero_lines=True, mu_sim=None)

• Call pyLLE.analyzedisp.AnalyzeDisp to get the dispersion of the resonator we want to simulate

LLEsolver.Setup(self)

• Set up the simulation for the Julia back-end. Save the two main dictionaries self.sim and self.res into a readable hdf5 file for Julia in the temporary location defined by the operating system.

LLEsolver.SolveTemporal(self, tol=0.001, maxiter=6,

step_factor=0.1)

• Call Julia to solve the LLE

LLEsolver.SolveSteadyState(self)

• Newton–Raphson method to find the root of the steady state equation

LLEsolver.RetrieveData(self)

• Load the output hdf5 saved by Julia and transform it into a user-friendly dictionary to be more suited for use with Python

LLEsolver.PlotCombPower(self,

do_matplotlib=False)

• Plot a figure with 3 subplots.

• •

Top subplot = map of the spectra for the steps taken by the LLE (step sub-sampled to be 1000)

• •

Middle subplot = temporal map of the intensity inside the resonator for the steps of the LLE

• •

Bottom subplot = normalized comb power

• •

Output

• –

f <obj>: matplotlib/plotly figure handle

• –

ax <obj>: matplotlib axes handle (if plotly, only returns f)

LLEsolver.PlotCombSpectra(self, ind,

f=None, ax=None, label=None, pwr=‘both’, do_matplotlib=False, plot=True)

• Plot the spectrum for a given index in the 1000 sub-sampled LLE steps

• •

Input

• –

ind <ind>: index in the LLE step for which the spectrum will be plotted

• –

f <obj>: matplotlib figure handle (if None, new figure)

• –

ax <obj>: matplotlib axes handle

• –

label <str>: label for the legend

• –

pwr <str>: ‘both’, ‘ring’, ‘wg’ depending on the spectra wanted (inside the ring, the waveguide, or both)

• •

Output

• –

freq <numpy.array>: frequency in Hz

• –

Sout <numpy.array>: spectral density of power in the waveguide (dBm)

• –

Sring <numpy.array>: spectral density of power in the ring (dBm)

• –

f <obj>: matplotlib/plotly figure handle

• –

ax <obj>: matplotlib axes handle (if plotly, only returns f)

LLEsolver.PlotSolitonTime(self, ind, f=None, ax=None,

label=None, do_matplotlib=False)

• Plot the spectra for a given index in the 1000 sub-sampled LLE step

• •

Input

• –

ind <ind>: index in the LLE step to plot the spectra

• –

f <obj>: matplotlib figure handle (if None, new figure)

• –

ax <obj>: matplotlib axe handle

• –

label <str>: label for the legend

• •

Output

• –

τ <obj>: Time in the resonator

• –

U <numpy.array>: Temporal electric field for the given step of the LLE

• –

f <obj>: matplotlib/plotly figure handle

• –

ax <obj>: matplotlib axes handle (if plotly, only returns f)

LLEsolver.SaveResults(self, fname,

path=‘./’)

Save the whole class using the pickle package, which makes it easy to load back the saved results

• •

Input

• –

fname <str>: name to save. The ‘.pkl’ extension will be added

• –

path <str>: path to save the results (defaults ‘./’)

5.2 _analyzedisp

5.2.1 AnalyzeDisp(self, **kwargs)

Call to analyze the dispersion of a simulated resonator. Initialization input. Everything is in SI ([]=optional):

• •

Input:

• –

f_center <float>: pump frequency

• –

file <str>: .txt file to load

• –

R <float>: radius of the resonator

• –

rM_fit <list>.: lower and upper bounds of mode to fit the dispersion

• –

rM_sim <list>.: lower and upper bounds of the modes to extrapolate for the simulation

• –

f <obj>: matplotlib figure handle

• –

ax <obj>: matplotlib axes handle

• –

label <list>: list of the string for each plot to be labeled

• –

plottype <str>: define the type of plot ‘all’ [defaults], ‘sim’, ‘fit’, ‘fem’, ‘ind’

AnalyzeDisp.GetDint(self)

Retrieve the dispersion of a resonator based on the frequency of resonance and azimuthal mode order. The data are fit using a cubic spline method

• •

Output:

• –

self.PrM_fit: scipy.interpolate object which fitted the data

• –

self.Dint_fit: fitted integrated dispersion for the simulated mode

• –

self.neff_pmp: effective index at the pump frequency

• –

self.ng_pmp: group index at the pump frequency

6 Acknowledgments

The authors acknowledge funding support from the DARPA DODOS and ACES programs, the NIST-on-a-chip program, and the Cooperative Research Agreement between the University of Maryland and NIST-CNST (award no. 70NANB10H193). They also thank Dr. Su-Peng Yu and Dr. Jared H. Strait for their helpful comments on the manuscript.