pyLLE: A Fast and User Friendly Lugiato-Lefever Equation Solver. (2024)

1. Introduction

The Lugiato-Lefever Equation (LLE), first developed to provide adescription of spatial dissipative structures in optical systems [1],has recently made a significant impact in the integrated photonicscommunity, where it has been adopted to help understand and predictKerr- mediated nonlinear optical phenomena such as parametric frequencycomb generation [2, 3] inside microresonators [4-6]. The LLE isessentially an application of the nonlinear Schrodinger equation (NLSE)[7] to a damped, driven Kerr nonlinear resonator, so that a periodicboundary condition is applied. Importantly, a slow-varying time envelopeis stipulated, resulting in a mean-field solution in which the fielddoes not vary within a round trip. This constraint, which differentiatesthe LLE from the more general Ikeda map [8], significantly simplifiescalculations while still providing excellent physical representation fora wide variety of systems. In particular, simulations based on the LLEformalism have enabled modeling that quantitatively agrees with reportedexperimental results on microcomb generation [9, 10] (e.g., in terms ofspectral bandwidth), and have also been central to theoretical studiesthat have provided better insight into novel nonlinear dynamics that canbe supported by Kerr nonlinear microresonators [11-13].

The great potential of microresonator frequency combs (microcombs)in a wide variety of applications [3, 14-22] suggests the need forefficient and widely accessible computational tools to more rapidlyfurther their development. Although LLE simulations are commonlyperformed by research groups working in the field, to our knowledge nofree software package for solving this equation in an easy and fast wayis currently available. Here, we introduce pyLLE, an open-source LLEsolver for microcomb modeling. It combines the user-friendliness of thePython programming language and the computational power of the Juliaprogramming language.

2. Software Specifications

NIST Operating Unit Physical Measurement Laboratory, Microsystems and Nanotechnology DivisionCategory Lugiato Lefever Equation SolverTargeted Users Scientist / Researcher / Engineer in Photonics/MetrologyOperating Systems Windows/Linux/MacOSProgramming Language Python>3.4, Julia 0.6.4Inputs/Outputs Input text file with resonator dispersion; output electric field dataDocumentation https://usnistgov.github.io/pyLLE/Accessibility N/ADisclaimer https://www.nist.gov/director/licensing

3. Lugiato-Lefever Equation Model

The Lugiato-Lefever equation is derived from the nonlinearSchrodinger equation, which can be written, if considering only thethird order nonlinearity and in the reference frame moving with thegroup velocity, as:

[mathematical expression not reproducible] (1)

with E (z, [tau]) being the complex electric field propagatingalong z during the time [tau], [alpha] the loss per unit length, [gamma]= [n.sub.2][[omega].sub.0]/([c.sub.0][A.sub.eff]) the effectivenonlinear coefficient, [n.sub.2] is the Kerr nonlinear coefficient(i.e., the nonlinear refractive index), [[omega].sub.0] the angularfrequency of nearest mode close to the pump frequency [[omega].sub.pmp],[A.sub.eff] the effective mode area in the resonator at the pumpfrequency, [[beta].sub.k] is the kth-order Taylor expansion coefficientof the dispersion [7], and [c.sub.0] the speed of light in vacuum.

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

[mathematical expression not reproducible] (2)

[mathematical expression not reproducible] (3)

where [E.sup.(m)] represents the field in the resonator for the mthround trip, [[phi].sub.0] the phase shift accumulated during a singleround trip, L = 2[pi]R the length of the resonator and [theta]theresonator-waveguide coupling coefficient.

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

[mathematical expression not reproducible] (4)

with t being the slow time and [tau] the fast time (i.e-[t.sub.R]/2 [less than or equal to] [tau] [less than or equal to][t.sub.R]/2), [t.sub.R] the round trip time, [alpha]' = [alpha]L +[theta] are the total round trip losses in the cavity and[[delta].sub.0] = [[omega].sub.pmp] - [[omega].sub.0] the detuning ofthe pump field with respect to the closest cold-cavity resonance[[omega].sub.0].

For systems with a complex dispersion profile, such as microringresonators, instead of using a Taylor expansion of the dispersion thatresults in a set of different [beta] coefficients (i.e., dispersion todifferent orders), it is often more convenient and accurate to work withthe integrated dispersion [D.sub.int] = [[omega].sub.[mu]] -([[omega].sub.0] + [D.sub.1][mu]), with [mu] an integer representing themode number relative to the pump resonance (for which [mu] = 0),[[omega].sub.[mu]] is the angular frequency of the [mu]th mode, and[D.sub.1] is the angular free-spectral range evaluated at the pumpedresonance (i.e., [D.sub.1] /2[pi] is the free-spectral range in units ofHz). This implies that one has to work in the frequency domain.Considering the approximation of the slowly-varying envelope, one canconsider the Fourier transform, FT, and inverse Fourier transform,[FT.sup.-1], relative to the fast time [tau]. Hence, Eq. (4) becomes theone implemented in the package presented here:

[mathematical expression not reproducible] (5)

This approach is similar to the one found in Hansson and Wabnitz[8], for simulating the Kerr dynamics across a large number of modesusing the Coupled Mode Theory (CMT) formalism. It is convenient indeedto work in the frequency domain, as it avoids the approximation oftruncating the dispersion beyond a certain order, and the resonancefrequencies can be found accurately using a Maxwell's equationeigenvalue solver. Moreover, working in the frequency domain allows theuse of well- developed computational libraries for executing the Fouriertransform operation in one iteration of the split-step Fouriertechnique, 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-moderesonator have been studied, including Coupled Mode Theory (CMT) and theLLE. However, the LLE provides a temporal approach which is well suitedto study Kerr non-linear resonators, especially in the context ofsoliton microcombs. Our fundamental goal in the pyLLE package is toprovide a fast, efficient, platform-independent and freely availabletool for simulating the LLE, to enable researchers with interest inmicrocombs to simulate their behavior without requiring in-depthknowledge of the LLE formalism or its numerical implementation.

The bulk of the computation is done through the Julia back-end. TheJulia language provides an efficient computational interface with BLAS(basic linear algebra subprograms), enabling reuse of the same FastFourier Transform plan (i.e plan_fft) throughout several iterations of aloop, 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 samesimulation of a Si3N4 ring resonator (described in more detail in thedocumentation 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 Pythonversion, and Matlab-based code. (1) We find that the simulation, inwhich the field behavior across 244 modes of the resonator iscalculated, is four times faster using pyLLE than the pure Pythonapproach, and almost twice as fast as a Matlab solver (Table 1). We notethat Matlab has a proprietary license, compared to the open-source pyLLEpackage.

As mentioned above, pyLLE is not only fast but highlyuser-friendly, thanks to the Python language which is becoming anincreasingly valuable resource in the scientific toolbox. Below wedescribe a full example of how to use pyLLE, which consists ultimatelyof 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 isset up and executed in pyLLE. We assume the pyLLE package has beeninstalled on the user local machine following the instructions availablein Ref. [23].

4.2.1 Setting up pyLLE

Let's start by importing the package:

import pyLLEOne has to define the resonator parameters.res = {'R': 23e-6, # ring radius in meters 'Qi' :1e6, # Intrinsic Q factor 'Qc': 1e6, # Coupling Q factor 'y': 1.55, # Effective nonlinear coefficient at the pump frequency (units of [W.sup.-1] [m.sup.-1]) 'dispfile': 'TestDispersion.txt',}

The parameters in the code above are mostly self-explanatory andare chosen to fit the physical device investigated. The dispfileparameter should be a raw text file, without header, in a "commaseparated value" (csv) format consisting of the integer azimuthalmode orders and corresponding resonance frequencies (in Hz) as the firstand second column, respectively. Importantly, the parameter names arethe ones commonly used in the LLE formalism, hence some follow the Greekalphabet. To improve the usability of pyLLE, one can define a parametereither by its Greek symbol or through its Latin name (e.g one can useeither [gamma] or gamma).

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

import numpy as npsim = {'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 '[delta][omega]_init': 2e9*2*np.pi,# Initial detuning of the pump in rad/s '[delta][omega]_end': -8e9*2*np.pi,# End detuning of the pump in rad/s '[mu]_sim': [-74,170], # Azimuthal mode numbers to simulate in the LLE on the left and right side of the pump '[mu]_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 withpump power Pin in the input coupling bus waveguide of 150 mW, with alinear detuning ramp of the pump from init to [delta][omega]-endrelative to the pump mode angular frequency, i.e., the mode closest tothe defined pump frequency f-pmp. The simulation length Tscan is inunits of number of round trips, as is most convenient in the LLEformalism. Two parameter ranges for the simulation frequency range haveto be defined, [mu]_fit, which determines the fit window for the rawdata found in dispfile, and [mu]_sim which is the number of modessimulated in the LLE. [mu]_sim can differ from [mu]_fit (i.e., thesimulation can be performed over a truncated range or can be expandedover a larger range through extrapolation).

Once the resonator and simulation parameters are set, one caneasily 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 methodsolver.Analyze has to be called, resulting in the plot of the integrateddispersion [D.sub.int], and returning the handle of the figure (and axesif using matplotlib).

In Fig. (2), one extrapolates (orange curve - fit dispersion)outside of the spectral window containing the raw dispersion data (greencircles), to enable a simulation over a broader spectral range (bluedashed line). One has to be careful about using this extrapolationfeature, to make sure that the extrapolated dispersion is a reasonablephysical representation of the dispersion. For example, extrapolationcan result in ripples in the integrated dispersion, creatingzero-crossings that are artifacts that will influence the LLE simulationresults.

4.2.3 Temporal Simulation

The core function of pyLLE is to solve the full temporalLugiato-Lefever Equation as described in Eq. (5). The solver implementedin this package is based on a Julia core called from python. Tointerface both languages easily, a .hdf5 file is created in a temporarylocation 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 thecurrent simulation, making it convenient to assess the time a simulationwill 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. Thesolver.PlotCombPower method (Fig. 3) provides an overview of themicrocomb behavior, plotting the spectra, time envelope, and comb powerinside the resonator for the different steps of the LLE, which are sub-sampled for faster display.

To obtain a better idea of the microcomb behavior for a givenlaser-cavity mode detuning, we plot the spectrum (Fig. 4) and thetemporal profile (Fig. 5) of the electric field in the resonator at afixed LLE step:

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

4.2.4 Steady-state Simulation

Instead of solving the full temporal LLE in Eq. (5), we can insteadconsider solutions to the steady-state LLE, i.e., considering theright-hand side of Eq. (5) equal to zero. Steady- state solutions to theLLE indicate what types of states are supported by the system, but donot necessarily give an indication of how to access such states inpractice.

Newton's method can be used to find steady-state solutions andis implemented in the method solver.SolveSteadyState. Although it givesfast results, the accuracy of such a solver remains questionablecompared 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 sameway 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, andcalculating the dispersion profile. Here it is not necessary to use themethod solver.Setup, as the Newton's method places relatively lowcomputational demands and thus everything is solved with Python.

However, the main difference with respect to the temporal solver isthe introduction of a fixed detuning [lambda][omega] of the pumprelative to the closest mode, so that:

solver.sim ['[lambda][omega]'] = -5e9*2*np.pi # thisvalue corresponds to the detuning at the end of the soliton step in thetemporal simulation

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

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 withthe 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 waveguidecomprising the resonator.

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

* sim <dict>

--Tscan <float>: length of the simulation (in units of roundtrip)

--mufit <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 detuningbut keep doing the simulation

* debug <bool>: Save a trace in a log file in the workingdirectory 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 theresonator we want to simulate LLEsolver.Setup(self)

Set up the simulation for the Julia back-end. Save the two maindictionaries self.sim and self.res into a readable hdf5 file for Juliain 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 auser-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 theresonator 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 returnsf)

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 LLEsteps

* Input

--ind <ind>: index in the LLE step for which the spectrumwill 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, thewaveguide, or both)

* Output

--freq <numpy.array>: frequency in Hz

--Sout <numpy.array>: spectral density of power in thewaveguide (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 returnsf)

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

--[tau] <obj>: Time in the resonator

--U <numpy.array>: Temporal electric field for the given stepof the LLE

--f < obj >: matplotlib/plotly figure handle

--ax <obj>: matplotlib axes handle (if plotly, only returnsf)

LLEsolver.SaveResults(self, fname, path='./') Save thewhole class using the pickle package, which makes it easy to load backthe saved results

* Input

--fname <str>: name to save. The '.pkl' extensionwill 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 thedispersion

--rMsim <list>.: lower and upper bounds of the modes toextrapolate for the simulation

--f < obj >: matplotlib figure handle

--ax < obj > : matplotlib axes handle

--label <list>: list of the string for each plot to belabeled

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

AnalyzeDisp.GetDint(self) Retrieve the dispersion of a resonatorbased on the frequency of resonance and azimuthal mode order. The dataare fit using a cubic spline method

* Output:

--self.PrMfit: scipy.interpolate object which fitted the data

--self.Dintfit: 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

Key words: Julia; micro-combs, non-linear optics; python.

Acknowledgments

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

6. References

[1] Lugiato LA, Lefever R (1987) Spatial Dissipative Structures inPassive Optical Systems. Physical Review Letters 58(21):2209--2211.https://doi.org/10.1103/PhysRevLett.58.2209

[2] Del'Haye P, Schliesser A, Arcizet O, Wilken T, HolzwarthR, Kippenberg TJ (2007) Optical frequency comb generation from amonolithic microresonator. Nature 450(7173):1214-1217.https://doi.org/10.1038/nature06401

[3] Kippenberg TJ, Gaeta AL, Lipson M, Gorodetsky ML (2018)Dissipative kerr solitons in optical microresonators. Science361(6402):eaan8083. https://doi.org/10.1126/science.aan8083

[4] Coen S, Randle HG, Sylvestre T, Erkintalo M (2012) Modeling ofoctave- spanning Kerr frequency combs using a generalized mean-fieldLugiato-Lefever model. Optics Letters 38(1):37-39.https://doi.org/10.1364/ol.38.000037

[5] Chembo YK, Menyuk CR (2013) Spatiotemporal Lugiato-Lefeverformalism for Kerr-comb generation in whispering-gallery-moderesonators. Physical Review A 87(5):053852-053854.https://doi.org/10.1103/PhysRevA.87.053852

[6] Lamont MRE, Okawachi Y, Gaeta AL (2013) Route to stabilizedultrabroadband microresonator-based frequency combs. Optics Letters38(18):3478-3481. https://doi.org/10.1364/OL.38.003478

[7] Agrawal GP (2007) Nonlinear Fiber Optics (Academic Press) 4thEd. https://doi.org/10.1016/b978-0-12-369516-1.x5000-6

[8] Hansson T, Wabnitz S (2015) Frequency comb generation beyondthe Lugiato- Lefever equation: multi-stability and super cavitysolitons. Journal of the Optical Society of America B 32(7):1259-1266.https://doi.org/10.1364/JOSAB.32.001259

[9] Xue X, Xuan Y, Liu Y, Wang PH, Chen S, Wang J, ELeaird D, Qi M,Weiner A (2015) Mode-locked dark pulse Kerr combs in normal-dispersionmicroresonators. Nature Photonics 9:594-600.https://doi.org/10.1038/nphoton.2015.137

[10] Li Q, Briles TC, Westly DA, Drake TE, Stone JR, Ilic BR,Diddams SA, Papp SB, Srinivasan K (2017) Stably accessingoctave-spanning microresonator frequency combs in the soliton regime.Optica 4(2):193-203. https://doi.org/10.1364/0PTICAA000193

[11] Kozyreff G (2012) Localized Turing patterns in nonlinearoptical cavities. Physica D: Nonlinear Phenomena 241(10):939-946.https://doi.org/10.1016Zj.physd.2012.02.007

[12] Luo R, Liang H, Lin Q (2016) Multicolor cavity soliton. OpticsExpress 24(15):16777-16787. https://doi.org/10.1364/0E.24.016777

[13] Moille G, Li Q, Kim S, Westly D, Srinivasan K (2018)Phased-locked two- color single soliton microcombs indispersion-engineered Si3N4 resonators. Optics Letters 43(12):2772-2775.https://doi.org/10.1364/0L.43.002772

[14] Papp SB, Beha K, Del'Haye P, Quinlan F, Lee H, Vahala KJ,Diddams SA (2014) Microresonator frequency comb optical clock. Optica1(1):10. https://doi.org/10.1364/0PTICA.L000010

[15] Suh MG, Yang QF, Yang KY, Yi X, Vahala KJ (2016)Microresonator soliton dual-comb spectroscopy. Science354(6312):600-603. https://doi.org/10.1126/science.aah6516

[16] Marin-Palomo P, Kemal JN, Karpov M, Kordts A, Pfeifle J,Pfeiffer MHP, Trocha P, Wolf S, Brasch V, Anderson MH, Rosenberger R,Vijayan K, Freude W, Kippenberg TJ, Koos C (2017) Microresonator- basedsolitons for massively parallel coherent optical communications. Nature546(7657):274-279. https://doi.org/10.1038/nature22387

[17] Spencer DT, Drake T, Briles TC, Stone J, Sinclair LC, FredrickC, Li Q, Westly D, Ilic BR, Bluestone A, Volet N, Komljenovic T, ChangL, Lee SH, Oh DY, Suh MG, Yang KY, Pfeiffer MHP, Kippenberg TJ, NorbergE, Theogarajan L, Vahala K, Newbury NR, Srinivasan K, Bowers JE, DiddamsSA, Papp SB (2018) An optical-frequency synthesizer using integratedphotonics. Nature 557(7703):81-85.https://doi.org/10.1038/s41586-018-0065-7

[18] Trocha P, Karpov M, Ganin D, Pfeiffer MHP, Kordts A, Wolf S,Krockenberger J, Marin-Palomo P, Weimann C, Randel S, Freude W,Kippenberg TJ, Koos C (2018) Ultrafast optical ranging usingmicroresonator soliton frequency combs. Science 359(6378):887-891.https://doi.org/10.1126/science.aao3924

[19] Suh MG, Vahala KJ (2018) Soliton microcomb range measurement.Science 359(6378):884-887. https://doi.org/10.1126/science.aao1968

[20] Dutt A, Joshi C, Ji X, Cardenas J, Okawachi Y, Luke K, GaetaAL, Lipson M (2018) On-chip dual-comb source for spectroscopy. ScienceAdvances 4(3):e1701858. https://doi.org/10.1126/sciadv.1701858

[21] Obrzud E, Rainer M, Harutyunyan A, Anderson MH, Liu J,Geiselmann M, Chazelas B, Kundermann S, Lecomte S, Cecconi M, Ghedina A,Molinari E, Pepe F, Wildi F, Bouchy F, Kippenberg TJ, Herr T (2019) Amicrophotonic astrocomb. Nature Photonics 13(1):31-35.https://doi.org/10.1038/s41566-018-0309-y

[22] Suh MG, Yi X, Lai YH, Leifer S, Grudinin IS, Vasisht G, MartinEC, Fitzgerald MP, Doppmann G, Wang J, Mawet D, Papp SB, Diddams SA,Beichman C, Vahala K (2019) Searching for exoplanets using amicroresonator astrocomb. Nature Photonics 13(1):25-30.https://doi.org/10.1038/s41566-018-0312-3

[23] pyLLE: a Fast and User Friendly Lugiato-Lefever EquationSolver, National Institute of Standards and Technology. Available athttps://github.com/usnistgov/pyLLE

Gregory Moille, Qing Li, and Xiyuan Lu are guest researchassociates in the Microsystems and Nanotechnology Division at NIST, andare postdoctoral scholars at the University of Maryland. KartikSrinivasan is a physicist in the Microsystems and NanotechnologyDivision at NIST and an adjunct professor at the University of Maryland.The National Institute of Standards and Technology is an agency of theU.S. Department of Commerce.

Gregory Moille (1,2), Qing Li (1,2,3), Xiyuan Lu (1,2), KartikSrinivasan (1,4)

(1) National Institute of Standards and Technology, Gaithersburg,MD 20899 USA

(2) Institute for Research in Electronics and Applied Physics andMaryland Nanocenter, University of Maryland, College Park, Maryland20742 USA

(3) Electrical and Computer Engineering, Carnegie MellonUniversity, Pittsburgh, PA 15213 USA

(4) Joint Quantum Institute, NIST/University of Maryland, CollegePark, Maryland 20742 USA

[emailprotected]

[emailprotected]

Software DOI: https://doi.org/10.18434/M32059

Software Version: 2.1

Accepted: April 29, 2019

Published: May 24, 2019

https://doi.org/10.6028/jres.124.012

(1) Certain commercial items are identified in this paper to fosterunderstanding. Such identification does not imply recommendation orendorsem*nt by NIST, nor does it imply that the items identified arenecessarily the best available for the purpose.

Caption: Fig. 1. Cartoon of a ring resonator coupled to a waveguidewhere a continuous wave pump is sent as an input and produces, through[[chi].sup.(3)]-mediated nonlinear interactions in the ring, a frequencycomb at the output. The coupling coefficient from the waveguide to thering and vice-versa is [theta] and the intrinsic resonator losses perunite of length is [alpha].

Caption: Fig. 2. Integrated dispersion retrieved through thedispfile file, with the raw values from the file (green circles), fitdispersion (orange solid line), and integrated dispersion that will beused in the solver (blue dashed line).

Caption: Fig. 3. (a) Comb spectra vs. LLE step. (b) Temporalenvelope of the electromagnetic field vs. LLE step. (c) Comb power(blue) vs. LLE step. The pump detuning during the simulation (orange) isalso plotted.

Caption: Fig. 4. Spectrum of the soliton comb obtained at the LLEstep = 570, inside the resonator (orange), and out-coupled into theaccess waveguide (blue).

Caption: Fig. 5. Electromagnetic envelope of the time-domainprofile of a single soliton state inside the resonator.

Caption: Fig. 6. Result of a steady-state simulation of the combspectrum within the resonator (orange) and out-coupled into the accesswaveguide (blue).

Table 1. Performance comparison of the same resonator LLEsimulation using three different software implementations.Matlab 2018a Pure Python pyLLE19 mins 45 min 11 mins