3. Basic Usage

3.1. A First Calculation

We’re now ready to start using NumBAT.

Let’s jump straight in and run a simple calculation. Later in the chapter, we go deeper into some of the details that we will encounter in this first example.

3.1.1. Tutorial 1 – Basic SBS Gain Calculation

Simulations with NumBAT are generally carried out using a python script file.

This example, contained in <NumBAT>tutorials/simo-tut_01-first_calc.py calculates the backward SBS gain for a rectangular silicon waveguide surrounded by air.

Move into the tutorials directory and then run the script by entering:

$ python3 simo-tut_01-first_calc.py

After a short while, you should see some values for the SBS gain printed to the screen. In many more tutorials in the subsequent chapters, we will meet much more convenient forms of output, but for now let’s focus on the steps involved in this basic calculation.

The sequence of operations (annotated in the source code below as Step 1, Step 2, etc) is:

1. Add the NumBAT install directory to Python’s module search path and then import

   the NumBAT python modules.

2. Set parameters to define the structure shape and dimensions.

3. Set parameters determining the range of electromagnetic and elastic modes to be solved.

4. Construct the waveguide with objects.Structure out of a number of materials.Material objects. The generated mesh is shown in the figure below.

5. Solve the electromagnetic problem at a given free space wavelength \(\lambda\). The function mode_calcs.calc_EM_modes() returns an object containing electromagnetic mode profiles, propagation constants, and potentially other data which can be accessed through various methods.

6. Display the propagation constants in units of \(\text{m}^{-1}\) of the EM modes using mode_calcs.kz_EM_all()

7. Obtain the effective index of the fundamental mode using mode_calcs.neff()

8. Identify the desired elastic wavenumber from the difference of the pump and Stokes propagation constants and solve the elastic problem. mode_calcs.calc_AC_modes() returns an object containing the elastic mode profiles, frequencies and potentially other data at the specified propagation constant k_AC.

9. Display the elastic frequencies in Hz using mode_calcs.nu_AC_all().

10. Calculate the total SBS gain, contributions from photoelasticity and moving boundary effects, and the elastic loss using integration.gain_and_qs(), and print them to the screen.

Note from this description that the eigenproblems for the electromagnetic and acoustic problems are framed in opposite senses. The electromagnetic problem finds the wavenumbers \(k_{z,n}(\omega)\) (or equivalently the effective indices) of the modes at a given free space wavelength (ie. at a specified frequency \(\omega=2\pi c/\lambda\)). The elastic solver, however, works in the opposite direction, finding the elastic modal frequencies \(\nu_n(q_0)\) at a given elastic propagation constant \(q_0\). While this might seem odd at first, it is actually the natural way to frame SBS calculations.

We emphasise again, that for convenience, the physical dimensions of waveguides are specified in nanometres. All other quantities in NumBAT are expressed in the standard SI base units.

Generated mesh for this example as displayed by guide.check_mesh().

Here’s the full source code for this tutorial:

print(nbapp.final_report())
""" Calculate the backward SBS gain for modes in a
    silicon waveguide surrounded in air.
"""

# Step 1
import sys
import numpy as np

sys.path.append("../backend/")

import numbat
import integration
import mode_calcs
import materials


# Naming conventions
# AC: acoustic
# EM: electromagnetic
# q_AC: acoustic wavevector

print('\n\nCommencing NumBAT tutorial 1')

# Step 2
# Geometric Parameters - all in nm.
lambda_nm = 1550  # Wavelength of EM wave in vacuum.
# Unit cell must be large to ensure fields are zero at boundary.
unitcell_x = 2.5*lambda_nm
unitcell_y = unitcell_x
# Waveguide widths.
inc_a_x = 300
inc_a_y = 280
# Shape of the waveguide.
inc_shape = 'rectangular'

# Step 3
# Number of electromagnetic modes to solve for.
num_modes_EM_pump = 20
num_modes_EM_Stokes = num_modes_EM_pump
# Number of acoustic modes to solve for.
num_modes_AC = 20
# The EM pump mode(s) for which to calculate interaction with AC modes.
# Can specify a mode number (zero has lowest propagation constant) or 'All'.
EM_ival_pump = 0
# The EM Stokes mode(s) for which to calculate interaction with AC modes.
EM_ival_Stokes = 0
# The AC mode(s) for which to calculate interaction with EM modes.
AC_ival = 'All'

# Step 4
# Use specified parameters to create a waveguide object.
# to save the geometry and mesh as png files in backend/fortran/msh/

nbapp = numbat.NumBATApp()


#wguide = nbapp.make_structure(unitcell_x, inc_a_x, unitcell_y, inc_a_y, inc_shape,
wguide = nbapp.make_structure(unitcell_x, inc_a_x, unitcell_y, inc_a_y, inc_shape,
                           material_bkg=materials.make_material("Vacuum"),
                           material_a=materials.make_material("Si_2016_Smith"),
                           lc_bkg=.1,  # in vacuum background
                           lc_refine_1=5.0,  # on cylinder surfaces
                           lc_refine_2=5.0)  # on cylinder center

# Note use of rough mesh for demonstration purposes by turning this line on.
# wguide.check_mesh()

# Explicitly remind ourselves what data we're using.
print('\nUsing material data: ', wguide.get_material('a'))

# Step 5
# Estimate expected effective index of fundamental guided mode.
n_eff = wguide.get_material('a').refindex_n-0.1

# Calculate the Electromagnetic modes of the pump field.
sim_EM_pump = wguide.calc_EM_modes(num_modes_EM_pump, lambda_nm, n_eff)

# Display the wavevectors of EM modes.
v_kz = sim_EM_pump.kz_EM_all()
print('\n k_z of electromagnetic modes [1/m]:')
for (i, kz) in enumerate(v_kz):
    print(f'{i:3d}  {np.real(kz):.4e}')

# Calculate the Electromagnetic modes of the Stokes field.
# For an idealised backward SBS simulation the Stokes modes are identical
# to the pump modes but travel in the opposite direction.
sim_EM_Stokes = mode_calcs.bkwd_Stokes_modes(sim_EM_pump)
# # Alt
# sim_EM_Stokes = wguide.calc_EM_modes(lambda_nm, num_modes_EM_Stokes, n_eff, Stokes=True)

# Step 6
# Find the EM effective index of the waveguide.
n_eff_sim = np.real(sim_EM_pump.neff(0))
print("\n Fundamental optical mode ")
print(" n_eff = ", np.round(n_eff_sim, 4))

# Acoustic wavevector
q_AC = np.real(sim_EM_pump.kz_EM(0) - sim_EM_Stokes.kz_EM(0))

print('\n Acoustic wavenumber (1/m) = ', np.round(q_AC, 4))

# Step 7
# Calculate Acoustic modes, using the mesh from the EM calculation.
sim_AC = wguide.calc_AC_modes(num_modes_AC, q_AC, EM_sim=sim_EM_pump)

# Print the frequencies of AC modes.
v_nu = sim_AC.nu_AC_all()
print('\n Freq of AC modes (GHz):')
for (i, nu) in enumerate(v_nu):
    print(f'{i:3d}  {np.real(nu)*1e-9:.5f}')

# Do not calculate the acoustic loss from our fields, instead set a Q factor.
set_q_factor = 1000.

# Step 8
# Calculate interaction integrals and SBS gain for PE and MB effects combined,
# as well as just for PE, and just for MB. Also calculate acoustic loss alpha.
SBS_gain_tot, SBS_gain_PE, SBS_gain_MB, linewidth_Hz, Q_factors, alpha = integration.gain_and_qs(
    sim_EM_pump, sim_EM_Stokes, sim_AC, q_AC, EM_ival_pump=EM_ival_pump,
    EM_ival_Stokes=EM_ival_Stokes, AC_ival=AC_ival, fixed_Q=set_q_factor)

# SBS_gain_tot, SBS_gain_PE, SBS_gain_MB are 3D arrays indexed by pump, Stokes and acoustic mode
# Extract those of interest as a 1D array:
SBS_gain_PE_ij = SBS_gain_PE[EM_ival_pump, EM_ival_Stokes, :]
SBS_gain_MB_ij = SBS_gain_MB[EM_ival_pump, EM_ival_Stokes, :]
SBS_gain_tot_ij = SBS_gain_tot[EM_ival_pump, EM_ival_Stokes, :]

# Print the Backward SBS gain of the AC modes.
print("\nContributions to SBS gain [1/(WM)]")
print("Acoustic Mode number | Photoelastic (PE) | Moving boundary(MB) | Total")

for (m, gpe, gmb, gt) in zip(range(num_modes_AC), SBS_gain_PE_ij, SBS_gain_MB_ij, SBS_gain_tot_ij):
    print(f'{m:8d}  {gpe:18.6e}  {gmb:18.6e}  {gt:18.6e}')


# Mask negligible gain values to improve clarity of print out.
threshold = 1e-3
masked_PE = np.where(np.abs(SBS_gain_PE_ij) > threshold, SBS_gain_PE_ij, 0)
masked_MB = np.where(np.abs(SBS_gain_MB_ij) > threshold, SBS_gain_MB_ij, 0)
masked_tot = np.where(np.abs(SBS_gain_tot_ij) > threshold, SBS_gain_tot_ij, 0)

print("\n Displaying gain results with negligible components masked out:")

print("AC Mode number | Photoelastic (PE) | Moving boundary(MB) | Total")
for (m, gpe, gmb, gt) in zip(range(num_modes_AC), masked_PE, masked_MB, masked_tot):
    print(f'{m:8d}  {gpe:18.6e}  {gmb:18.6e}  {gt:18.6e}')

print(nbapp.final_report())

In the next three chapters, we meet many more examples that show the different capabilities of NumBAT and provided comparisons against analytic and experimental results from the literature.

For the remainder of this chapter, we will explore some of the details involved in specifying a wide range of waveguide structures.

3.2. General Simulation Procedures

Simulations with NumBAT are generally carried out using a python script file. This file is kept in its own directory which may or may not be within your NumBAT tree. All results of the simulation are automatically created within this directory. This directory then serves as a complete record of the calculation. Often, we will also save the simulation objects within this directory for future inspection, manipulation, plotting, etc.

These files can be edited using your choice of text editor (for instance nano or vim) or an IDE (for instance MS Visual Code or pycharm) which allow you to run and debug code within the IDE.

To save the results from a simulation that are displayed upon execution (the print statements in your script) use:

$ python3 ./simo-tut_01-first_calc.py | tee log-simo.log

To have direct access to the simulation objects upon the completion of a script use:

$ python3 -i ./simo-tut_01-first_calc.py

This will execute the python script and then return you into an interactive python session within the terminal. This terminal session provides the user experience of an ipython type shell where the python environment and all the simulation objects are in the same state as in the script when it has finished executing. In this session you can access the docstrings of objects, classes and methods. For example:

>>> from pydoc import help
>>> help(objects.Struct)

where we have accessed the docstring of the Struct class from objects.py.

3.3. Script Structure

As with our first example above, most NumBAT scripts proceed with a standard structure:

• importing NumBAT modules

• defining materials

• defining waveguide geometries and associating them with material properties

• solving electromagnetic and acoustic modes

• calculating gain and other derived quantities

The following section provides some information about specifying material properties and waveguide structures, as well as the key parameters for controlling the finite-element meshing. Information on how to add new structures to NumBAT is provided in Making New Mesh.

3.4. Materials

In order to calculate the modes of a structure we must specify the acoustic and optical properties of all constituent materials.

In NumBAT, this data is read in from human-readable .json files, which are stored in the directory <NumBAT>/backend/material_data.

These files not only provide the numerical values for optical and acoustic variables, but provide links to the origin of the data. Often they are taken from the literature and the naming convention allows users to select from different parameter values chosen by different authors for the same nominal material.

The intention of this arrangement is to create a library of materials that can we hope can form a standard amongst the research community. They also allow users to check the sensitivity of their results on particular parameters for a given material.

At present, the library contains the following materials:

• Vacuum (or air)

  • Vacuum

• The chalcogenide glass Arsenic tri-sulfide

  • As2S3_2016_Smith

  • As2S3_2017_Morrison

  • As2S3_2021_Poulton

• Fused silica

  • SiO2_2013_Laude

  • SiO2_2015_Van_Laer

  • SiO2_2016_Smith

  • SiO2_2021_Smith

  • SiO2_smf28.json

  • SiO2GeO2_smf28.json

• Silicon

  • Si_2012_Rakich

  • Si_2013_Laude

  • Si_2015_Van_Laer

  • Si_2016_Smith

  • Si_2021_Poulton

  • Si_test_anisotropic

• Silicon nitride

  • Si3N4_2014_Wolff

  • Si3N4_2021_Steel

• Gallium arsenide

  • GaAs_2016_Smith

• Germanium

  • Ge_cubic_2014_Wolff

• Lithium niobate

  • LiNbO3_2021_Steel

  • LiNbO3aniso_2021_Steel

Materials can easily be added to this library by copying any of these files as a template and modifying the properties to suit. The Si_test_anisotropic file contains all the variables that NumBAT is setup to read. We ask that stable parameters (particularly those used for published results) be added to the NumBAT git repository using the same naming convention.

3.5. Waveguide Geometries

NumBAT encodes different waveguide structures through finite element meshes constructed using the .geo language used by the open source tool Gmsh. Most users will find they can construct all waveguides of interest using the existing templates. However, new templates can be added by adding a new .geo file to the <NumBAT>/backend/fortran/msh directory and making a new subclass of the UserGeometryBase class in the <NumBAT>/backend/msh/user_meshes.py file. Interested users can get in touch with <michael.steel@mq.edu.au>. All the builtin examples below are constructed in the same fashion in a parallel builtin_meshes.py file and can be used as models for your own designs.

The following figures give some examples of how material types and physical dimensions are represented in the mesh geometries. In particular, for each structure template, they identify the interpretation of the dimensional parameters (inc_a_x, slab_b_y, etc), material labels (material_a, material_b etc), and the grid refinement parameters (lc_bkg, lc_refine_1, lc_refine_2, etc). The captions for each structure also identify the mesh geometry template files in the directory <NumBAT>/backend/fortran/msh with filenames of the form <prefix>_msh_template.geo which define the structures and can give ideas for developing new structure files.

The NumBAT code for creating all these structures can be found in <NumBAT>/docs/source/images/make_meshfigs.py.

3.5.1. Single inclusion waveguides with surrounding medium

These structures consist of a single medium inclusion (mat_a) with a background material (mat_bkg). The dimensions are set with inc_a_x and inc_a_y.

Rectangular waveguide using shape rectangular (template oneincl_msh).

Elliptical waveguide using shape circular (template oneincl_msh).

Triangular waveguide using shape triangular.

3.5.2. Double inclusion waveguides with surrounding medium

These structures consist of a pair of waveguides with a single common background material. The dimesions are set by inc_a_x/inc_a_y and inc_b_x/inc_b_y. They are separated horintally by two_inc_sep and the right waveguide has a vertical offset of y_off.

Coupled rectangular waveguides using shape rectangular (template twoincl_msh).

Coupled circular waveguides using shape circular (template twoincl_msh). There appers to be a bug here!

3.5.3. Rib waveguides

These structures consist of a rib on one or more substrate layers with zero to two coating layers.

A conventional rib waveguide using shape rib (template rib).

A coated rib waveguide using shape rib_coated (template rib_coated).

A rib waveguide on two substrates using shape rib_double_coated (template rib_double_coated).

3.5.4. Engineered rib waveguides

These are examples of more complex rib geometries. These are good examples to study in order to make new designs using the user-specified waveguide and mesh mechanism.

A trapezoidal rib structure using shape trapezoidal_rib.

A supported pedestal structure using shape pedestal.

3.5.5. Slot waveguides

These slot waveguides can be used to enhance the horizontal component of the electric field in the low index region by the ‘slot’ effect.

A slot waveguide using shape slot (material_a is low index) (template slot).

A coated slot waveguide using shape slot_coated (material_a is low index) (template slot_coated).

3.5.6. Layered circular waveguides

These waveguides consist of a set of concentric circular rings of a desired number of layers in either a square or circular outer domain. Note that inc_a_x specifies the innermost diameter. The subsequent parameters inc_b_x, inc_c_x, etc specify the annular thickness of each succssive layer.

A many-layered concentric structure using shape onion (template onion).

A two-layered concentric structure with background using shape onion2 (template onion2).

A three-layered concentric structure with background using shape onion3 (template onion3).

A many-layered concentric structure with a circular outer boundary using shape circ_onion (template circ_onion).

A two-layered concentric structure with a circular outer boundary using shape circ_onion2 (template circ_onion2).

A three-layered concentric structure with a circular outer boundary using shape circ_onion3 (template circ_onion3).

3.6. User-defined waveguide geometries

Users may incorporate their own waveguide designs fully into NumBAT with the following steps. The triangular built-in structure is a helpful model to follow.

• Create a new gmsh template .geo file to be placed in <NumBAT>/backend/msh that specifies the general structure. Start by looking at the structure of triangular_msh_template.geo and some other files to get an idea of the general structure. We’ll suppose the file is called mywaveguide_msh_template.geo and the template name is thus mywaveguide.

  When designing your template, please ensure the following:

  • That you use appropriate-sized parameters for all significant dimensions. This makes it easier to determine if the template structure has the right general shape, even though the precise dimensions will usually be changed through NumBAT calls.

  • That all Line elements are unique. In other words do not create two Line objects joining the same two points. This will produce designs that look correct, but lead to poorly formed meshes that will fail when NumBAT runs.

  • That all Line Loop elements defining a particular region are defined with the same handedness. The natural choice is to go around the loop anti-clockwise. Remember to include a minus sign for any line element that is traversed in the backwards sense.

  • That all regions that define a single physical structure with a common material are grouped together as a single Surface and then Physical Surface.

  • That the outer boundary is grouped as a Line Loop and then a Physical Line.

  • That the origin of coordinates is placed in a sensible position, such as a symmetry point close to where you expect the fundamental mode fields to be concentrated. This doesn’t actually affect NumBAT calculations but will produce more natural axis scales in output plots.

  You can see all examples of these principles followed in the mesh structures supplied with NumBAT.

• If this is your first, user-defined geometry, copy the file ‘’user_waveguides.json_template`` in <NumBAT>/backend/msh/ to user_waveguides.json in the same directory. This will ensure that subsequent git pull commands will not overwrite your work.

• Open the file user_waveguides.json and add a new dictionary element for your new waveguide, copying the general format of the pre-defined example entries.

• Fill in values for the wg_impl (the name of the python file implementing your waveguide geometry), wg_class (the name of the python class corresponding to your waveguide) and inc_shape (the waveguide template name) fields.

  • The value of inc_shape will normally be the your chosen template name, in this case mywaveguide. The other parameters can be chosen as you wish. It is natural to choose a class name which matches your template name, so perhaps MyWaveguide. However, depending on the number of geometries you create, it may be convenient to store all your classes in one python file so the filename for wg_impl may be the same for all your entries.

  • The active field allows a waveguide to be disabled if it is not yet fully working and you wish to use other NumBAT models in the meantime. You must set active to True of 1 in order to test your waveguide model.

  • Then save and close this file.

• Open or create the python file you just specified in the wg_impl field. This file must be placed in the <NumBAT>/backend/msh directory.

  • The python file must include the import line from usermesh import UserGeometryBase.

  • Create your waveguide class MyWaveguide by subclassing the UserGeometryBase class and adding init_geometry and apply_parameters methods using the Triangular class in builtin_meshes.py as a model. Both methods must take only self as arguments.

  • The init_geometry method specifies a few values including the name of the template .geo file, the number of distinct waveguide components and a short description.

  • The apply_parameters method is the mechanism for associating standard NumBAT symbols like inc_a_x, slab_a_y, etc with actual dimensions in your .geo file. This is done by string substitution of unique expressions in your .geo file using float values evaluated from the NumBAT parameters. Again, look at the examples in the Triangular class to see how this works.

  • Optionally, you may also add a draw_mpl_frame method. This provides a mechanism to draw waveguide outlines onto mode profile images and will be called automatically any time an electromagnetic or elastic mode profile is generated. The built-in waveguides Circular, Rectangular and TwoIncl provide good models for this method.

Designing and implementing a few waveguide structure should not be a daunting task but some steps can be confusing the first time round. If you hit any hiccups or have suggestions for trouble-shooting, please let us know.

3.7. Mesh parameters

The parameters lc_bkg, lc_refine_1, lc_refine_2 labelled in the above figures control the fineness of the FEM mesh and are set when constructing the waveguide, as discussed in the next chapter. The first parameter lc_bkg sets the reference background mesh size, typically as a fraction of the length of the outer boundary edge. A larger lc_bkg yields a coarser mesh. Reasonable starting values are lc_bkg=0.1 (10 mesh points on the outer boundary) to lc_bkg=0.05 (20 mesh points on the outer boundary).

As well as setting the overall mesh scale with lc_bkg, one can also refine the mesh near interfaces and near select points in the domain, as may be observed in the figures in the previous section. This helps to increase the mesh resolution in regions where there the electromagnetic and acoustic fields are likely to be strong and/or rapidly varying. This is achieved using the lc_refine_n parameters as follows. At the interface between materials, the mesh is refined to have characteristic length lc_bkg/lc_refine_1, therefore a larger lc_refine_1 gives a finer mesh by a factor of lc_refine_1 at these interfaces. The meshing program Gmsh automatically adjusts the mesh size to smoothly transition from a point that has one mesh parameter to points that have other meshing parameters. The mesh is typically also refined in the vicinity of important regions, such as in the center of a waveguide, which is done with lc_refine_2, which analogously to lc_refine_1, refines the mesh size at these points as lc_bkg/lc_refine_2.

For more complicated structures, there are additional lc_refine_<n> parameters. To see their exact function, look for these expressions in the particular .geo file.

Choosing appropriate values of lc_bkg, lc_refine_1, lc_refine_2 is crucial for NumBAT to give accurate results. The appropriate values depend strongly on the type of structure being studied, and so we strongly recommended carrying out a convergence test before delving into new structures (see Tutorial 5 for an example) starting from similar parameters as used in the tutorial simulations.

As will as giving low accuracy, a structure with too coarse a mesh is often the cause of the eigensolver failing to converge in which case NumBAT will terminate with an error. If you encounter such an error, try the calculation again with a slightly smaller value for lc_bkg, or slightly higher values for the lc_refine_n parameters.

On the other hand, it is wise to begin with relatively coarse meshes. It will be apparent that the number of elements scales roughly quadratically with the lc_refine parameters and so the run-time increases rapidly as the mesh becomes finer. For each problem, some initial experimentation to identify a mesh resolution that gives reasonable convergence in acceptable simulation is usually worthwhile.

3.8. Viewing the mesh

When NumBAT constructs a waveguide, the template geo file is converted to a concrete instatiation with the lc_refine and geometric parameters adjusted to the requested values. This file is then converted into a gmsh .msh file. When exploring new structures and their convergence behaviour, it is a very good idea to view the generated mesh frequently.

You can examine the resolution of your mesh by calling the plot_mesh(<prefix>) or check_mesh() methods on a waveguide Structure object. The first of these functions saves a pair of images of the mesh to a <prefix>-mesh.png file in the local directory which can be viewed with your preferred image viewer; the second opens the mesh in a gmsh window (see Tutorial 1 above).

In addition, the .msh file generated by NumBAT in any calculation is stored in <NumBAT>/backend/fortran/msh/build and can be viewed by running the command

gmsh <msh_filename>.msh

In some error situations, NumBAT will explicitly suggest viewing the mesh and will print out the required command to do so.