4. Tutorial

4.1. Introduction

This chapter provides a sequence of graded tutorials for learning NumBAT, exploring its applications and validating it against literature results and analytic solutions where possible. Before attempting your own calculations with NumBAT, we strongly advise working through the sequence of tutorial exercises which are largely based on literature results.

You may then choose to explore relevant examples drawn from a recent tutorial paper by Dr Mike Smith and colleagues, and a range of other literature studies, which are provided in the following two chapters, JOSA-B Tutorial Paper and Literature Examples.

4.2. Some Key Symbols

As far as practical we use consisent notation and symbols in the tutorial files. The following list introduces a few commonly encountered ones. Note that with the exception of the free-space wavelength \(\lambda\) and the spatial dimensions of waveguide structures, which are both specified in nanometres (nm), all quantities in NumBAT should be expressed in the standard SI units. For example, elastic freqencies \(\nu\) are expressed in Hz, not GHz.

lambda_nm

This is the free-space optical wavelength \(\lambda\) satisfying \(\lambda = 2\pi c/\omega\), where \(c\) is the speed of light and \(\omega\) is the angular frequency. For convenience, this parameter is specified in nm.

For most examples, we use the conventional value \(\lambda\) =1550 nm.

omega, omega_EM, om_EM

This is the electromagnetic angular frequency \(\omega = 2 \pi c/\lambda\) specified in \(\mathrm{rad.s}^{-1}\).

k, beta, k_EM

This is the electromagnetic wavenumber or propagation constant \(k\) or \(\beta\), specified in \(\mathrm{m}^{-1}\).

neff, n_eff

This is the electromagnetic modal effective index \(\bar{n}=ck/\omega\), which is dimensionless.

nu, nu_AC

This is the acoustic frequency \(\nu\) specified in Hz.

Omega, Omega_AC, Om_AC

This is the acoustic angular frequency \(\Omega = 2 \pi \nu\) specified in \(\mathrm{rad.s}^{-1}\).

q, q_AC

This is the acoustic wavenumber or propagation constant \(q=v_{ac} \Omega\), where \(v_{ac}\) is the phase speed of the wave. The acoustic wavenumber is specified in \(\mathrm{m}^{-1}\).

m

This is an integer corresponding to the mode number \(m\) of an electromagnetic mode \(\vec E_m(\vec r)\) or an acoustic mode \(\vec u_m(\vec r)\).

For both electromagnetic and acoustic modes, counting of modes begins with m=0.

For the electromagnetic problem in which frequency/free-space wavelength is the independent variable, the \(m=0\) mode has the highest effective index \(\bar{n}\) and highest wavenumber \(k\) of any mode for a given angular frequency \(\omega\).

For the acoustic problem, the wavenumber \(q\) is the independent variable and we solve for frequency \(\nu=\Omega/(2\pi)\). The \(m=0\) mode has the lowest frequency \(\nu\) of any mode for a given wavenumber \(q\).

inc_a_x, inc_a_y, inc_b_x, inc_b_y, slab_a_x, slab_a_y, … etc

These are dimensional parameters specifying the lengths of different aspects of a given structure: rib height, fibre radius etc. For convenience, these parameters are specified in nm.

4.3. Tutorials

We now walk through a number of simple simulations that demonstrate the basic use of NumBAT located in the <NumBAT>/tutorials directory.

We will meet a significant number of NumBAT functions in these tutorials, though certainly not all. The full Python interface is documented in the section _chap-pythonapi-label.

4.3.1. Tutorial 2 – SBS Gain Spectra

The first example we met in the previous chapter only printed numerical data to the screen with no graphical output. This example, contained in <NUMBAT>/tutorials/simo-tut_02-gain_spectra-npsave.py considers the same silicon-in-air structure but adds plotting of fields, gain spectra and techniques for saving and reusing data from earlier calculations.

As before, move to the <NUMBAT>/tutorials directory, and then run the calculation by entering:

$ python3 simo-tut_02-gain_spectra-npsave.py

Or you can take advantage of the Makefile provided in the directory and just type:

$ make tut02

Some of the tutorial problems can take a little while to run, especially if your computer is not especially fast. To save time, you can run most problems with a coarser mesh at the cost of somewhat reduced accuracy, by adding the flag fast=1 to the command line:

$ python3 simo-tut_02-gain_spectra-npsave.py fast=1

Or using the makefile technique, simply

$ make ftut02

The calculation should complete in a minute or so. You will find a number of new files in the tutorials directory beginning with the prefix tut_02 (or ftut_02 if you ran in fast mode).

4.3.1.1. Gain Spectra

The Brillouin gain spectra and plotted using the functions integration.gain_and_qs() and plotting.plot_gain_spectra(). The results are contained in the file tut_02-gain_spectra.png which can be viewed in any image viewer. On Linux, for instance you can use

$ eog tut_02_gain_spectra.png

to see this image:

Gain spectrum in tut_02-gain_spectra.png showing gain due to the photoelastic effect, gain due to moving boundary effect, and the total gain. The numbers near the main peaks identify the acoustic mode associated with the resonance.

Note how the different contributions from the photoelastic and moving-boundary effects are visible. In some cases, the total gain (blue) may be less than one or both of the separate effects if the two components act with opposite sign. JOSA-B Tutorial Paper and Literature Examples. (See Literature example 1 in the chapter Literature Examples for an interesting cexample of this phenomenon.)

Note also that prominent resonance peaks in the gain spectrum are labelled with the mode number \(m\) of the associated acoustic mode. This makes it easy to find the spatial profile of the most relevant modes (see below).

4.3.1.2. Mode Profiles

The choice of parameters for plot_gain_spectra() has caused several other files to be generated showing a zoomed-in version near the main peak, and the whole spectrum on \(\log\) and dB scales:

Zoom-in of the gain spectrum in the previous figure in the file tut_02-gain_spectra_zoom.png .

Gain spectrum viewed on a log scale in the field tut_02-gain_spectra-logy.png .

This example has also generated plots of some of the electromagnetic and acoustic modes that were found in solving the eigenproblems. These are created using the calls to plotting.plot_mode_fields() and stored in the sub-directory tut_02-fields.

Note that a number of useful parameters are also displayed at the top-left of each mode profile. These parameters can also be extracted using a range of function calls on a Mode object (see the API docs).

Electric field profile of the fundamental (\(m=0\)) optical mode profile stored in tut_02-fields/EM_E_field_00.png. The figures shows the modulus of the whole electric field \(|{\vec E}|^2\), a vector plot of the transverse field \({\vec E}_t=(E_x,E_y)\), and the three components of the electric field. NumBAT chooses the phase of the mode profile such that the transverse components are real. Note that the \(E_z\) component is \(\pi/2\) out of phase with the transverse components. (Since the structure is lossless, the imaginary parts of the transverse field, and the real part of \(E_z\) are zero).

Magnetic field profile of the fundamental (\(m=0\)) optical mode profile showing modulus of the whole magnetic field \(|{\vec H}|^2\), vector plot of the transverse field \({\vec H}_t=(H_x,H_y)\), and the three components of the magnetic field. Note that the \(H_z\) component is \(\pi/2\) out of phase with the transverse components.

Displacement field \(\vec u(\vec r)\) of the \(m=3\) acoustic mode with gain dominated by the moving boundary effect (green curve in gain spectra). Note that the frequency of \(\Omega/(2\pi)=11.99\) GHz (listed in the upper-left corner) corresponds to the first peak in the gain spectrum.

Displacement field \(\vec u(\vec r)\) of the \(m=4\) acoustic mode with gain dominated by the photo-elastic effect (red curve in gain spectra). Note that the frequency of \(\Omega/(2\pi)=13.45\) GHz corresponds to the second peak in the gain spectrum.

4.3.1.3. Miscellaneous comments

Here are some further elements to note about this example:

1. When using the fast= mode, the output data and fields directory begin with ftut_02 rather than tut_02.

2. It is frequently useful to be able to save and load the results of simulations to adjust plots without having to repeat the entire calculation. Here the flag recalc_fields determines whether the calculation should be done afresh and use previously saved data. This is performed using the save_simulation() and load_simulation() calls.

3. Plots of the modal field profiles are obtained using plotting.plot_mode_fields. Both electric and magnetic fields can be selected using EM_E or EM_H as the value of the EM_AC argument. The mode numbers to be plotted is specified by ivals. These fields are stored in a folder tut_02-fields/ within the tutorial folder. Later we will see how an alternative approach in which we extract a Mode object from a Simulation which represent a single mode that is able to plot itself. This can be more convenient.

4. The overall amplitude of the modal fields is arbitrary. In NumBAT, the maximum value of the electric field is set to be 1.0, and this may be interpreted as a quantity in units of V/m, \(\sqrt{\mathrm{W}}\) or other units as desired. Importantly, the plotted magnetic field \(\vec H(\vec r)\) is multiplied by the impedance of free space \(Z_0=\sqrt{\mu_0/\epsilon_0}\) so that \(Z_0 \vec H(\vec r)\) and \(\vec E(\vec r)\) have the same units, and the relative amplitudes between the electric and magnetic field plots are meaningful.

5. The suppress_imimre option suppresses plotting of the \(\text{Im}[x]\), \(\text{Im}[y]\) and \(\text{Re}[z]\) components of the fields which in a lossless non-leaky problem should normally be zero at all points and therefore not useful to plot.

6. By default, plots are exported as png format. Pass the option pdf_png=pdf to plot functions to generate a pdf output.

7. Plots of both spectra and modes are generated with a best attempt at font sizes, line widths etc, but the range of potential cases make it impossible to find a selection that works in all cases. Most plot functions therefore support the passing of a plotting.Decorator object that can vary the settings of these parameters and also pass additional commands to write on the plot axes. See the plotting API for details. This should be regarded as a relatively advanced NumBAT feature.

8. Vector field plots often require tweaking to get an attractive set of vector arrows. The quiver_points option controls the number of arrows drawn along each direction.

9. The plot functions and the Decorator class support many options. Consult the API chapter for details on how to fine tune your plots.

The full code for this simulation is as follows:

print(nbapp.final_report())
"""
    NumBAT Tutorial 2

    Calculate the backward SBS gain spectra of a silicon waveguide surrounded in air.

    Show how to save simulation objects (eg. EM mode calcs) to expedite the process
    of altering later parts of simulations.

    Show how to implement integrals in python and how to load data from Comsol.
"""

import sys
import numpy as np

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


import numbat
import plotting
import integration
import mode_calcs
import materials



# Geometric Parameters - all in nm.
lambda_nm = 1550
unitcell_x = 2*lambda_nm
unitcell_y = unitcell_x
inc_a_x = 300
inc_a_y = 280
inc_shape = 'rectangular'

num_modes_EM_pump = 20
num_modes_EM_Stokes = num_modes_EM_pump
num_modes_AC = 25
EM_ival_pump = 0
EM_ival_Stokes = 0
AC_ival = 'All'

# choose between faster or more accurate calculation
if len(sys.argv) > 1 and sys.argv[1] == 'fast=1':
    prefix = 'ftut_02'
    refine_fac = 1
else:
    prefix = 'tut_02'
    refine_fac = 5

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

nbapp = numbat.NumBATApp(prefix)

# Use of a more refined mesh to produce field plots.
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, lc_refine_1=5.0*refine_fac, lc_refine_2=5.0*refine_fac)


# wguide.check_mesh()

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

recalc_fields = True     # run the calculation from scratch
# recalc_fields=False   # reuse saved fields from previous calculation

if recalc_fields:
    # Calculate Electromagnetic modes.
    sim_EM_pump = wguide.calc_EM_modes(num_modes_EM_pump, lambda_nm, n_eff)
    sim_EM_Stokes = mode_calcs.bkwd_Stokes_modes(sim_EM_pump)

    print('\nSaving EM fields')
    sim_EM_pump.save_simulation('tut02_wguide_data')
    sim_EM_Stokes.save_simulation('tut02_wguide_data2')
else:
    # Once npz files have been saved from one simulation run,
    # set recalc_fields=True to use the saved data
    sim_EM_pump = mode_calcs.load_simulation('tot02_wguide_data')
    sim_EM_Stokes = mode_calcs.load_simulation('tot02_wguide_data2')

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

# Plot the E fields of the EM modes fields - specified with EM_AC='EM_E'.
# Zoom in on the central region (of big unitcell) with xlim_, ylim_ args,
# which specify the fraction of the axis to remove from the plot.
# For instance xlim_min=0.4 will remove 40% of the x axis from the left outer edge
# to the center. xlim_max=0.4 will remove 40% from the right outer edge towards the center.
# This leaves just the inner 20% of the unit cell displayed in the plot.
# The ylim variables perform the equivalent actions on the y axis.

# Let's plot fields for only the fundamental (ival = 0) mode.

print('\nPlotting EM fields')
# Plot the E field of the pump mode
plotting.plot_mode_fields(sim_EM_pump, xlim_min=0.4, xlim_max=0.4, ylim_min=0.4,
                          ylim_max=0.4, ivals=[EM_ival_pump], contours=True,
                          EM_AC='EM_E', ticks=True)

# Plot the H field of the pump mode
plotting.plot_mode_fields(sim_EM_pump, xlim_min=0.4, xlim_max=0.4, ylim_min=0.4,
                          ylim_max=0.4, ivals=[EM_ival_pump], contours=True,
                          EM_AC='EM_H', ticks=True)

# Calculate the EM effective index of the waveguide.
n_eff_sim = np.real(sim_EM_pump.neff(0))
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))

if recalc_fields:
    # Calculate and save acoustic modes.
    sim_AC = wguide.calc_AC_modes(num_modes_AC, q_AC, EM_sim=sim_EM_pump)

    print('Saving AC fields')
    sim_AC.save_simulation('tot02_wguide_data_AC')
else:
    sim_AC = mode_calcs.load_simulation('tot02_wguide_data_AC')

# 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}')

# The AC modes are calculated on a subset of the full unitcell,
# which excludes vacuum regions, so there is usually no need to restrict the area plotted
# with xlim_min, xlim_max etc.

print('\nPlotting acoustic modes')
plotting.plot_mode_fields(sim_AC, contours=True,
                          ticks=True, quiver_points=20, ivals=range(10))

if recalc_fields:
    # Calculate the acoustic loss from our fields.
    # Calculate interaction integrals and SBS gain for PE and MB effects combined,
    # as well as just for PE, and just for MB.
    SBS_gain, 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)
    # Save the gain calculation results
    np.savez('tut02_wguide_data_AC_gain', SBS_gain=SBS_gain, SBS_gain_PE=SBS_gain_PE,
             SBS_gain_MB=SBS_gain_MB, linewidth_Hz=linewidth_Hz)
else:
    npzfile = np.load('wguide_data_AC_gain.npz', allow_pickle=True)
    SBS_gain = npzfile['SBS_gain']
    SBS_gain_PE = npzfile['SBS_gain_PE']
    SBS_gain_MB = npzfile['SBS_gain_MB']
    linewidth_Hz = npzfile['linewidth_Hz']

# The following function shows how integrals can be implemented purely in python,
# which may be of interest to users wanting to calculate expressions not currently
# included in NumBAT. Note that the Fortran routines are much faster!
# Also shows how field data can be imported (in this case from Comsol) and used.
comsol_ivals = 5  # Number of modes contained in data file.
SBS_gain_PE_py, alpha_py, SBS_gain_PE_comsol, alpha_comsol = integration.gain_python(
    sim_EM_pump, sim_EM_Stokes, sim_AC, q_AC, 'comsol_ac_modes_1-5.dat',
    comsol_ivals=comsol_ivals)

# Print the PE contribution to gain SBS gain of the AC modes.
print("\n Displaying results of first five modes with negligible components masked out")
# Mask negligible gain values to improve clarity of print out.
threshold = -1e-3
masked_PE = np.ma.masked_inside(
    SBS_gain_PE[EM_ival_pump, EM_ival_Stokes, :comsol_ivals], 0, threshold)
print("SBS_gain [1/(Wm)] PE NumBAT default (Fortran)\n", masked_PE)
masked = np.ma.masked_inside(
    SBS_gain_PE_py[EM_ival_pump, EM_ival_Stokes, :], 0, threshold)
print("SBS_gain [1/(Wm)] python integration routines \n", masked)
masked = np.ma.masked_inside(
    SBS_gain_PE_comsol[EM_ival_pump, EM_ival_Stokes, :], 0, threshold)
print("SBS_gain [1/(Wm)] from loaded Comsol data \n", masked)

# Construct the SBS gain spectrum, built from Lorentzian peaks of the individual modes.
freq_min = np.real(sim_AC.nu_AC_all()[0]) - 2e9  # Hz
freq_max = np.real(sim_AC.nu_AC_all()[-1]) + 2e9  # Hz
plotting.plot_gain_spectra(sim_AC, SBS_gain, SBS_gain_PE, SBS_gain_MB, linewidth_Hz,
                           EM_ival_pump, EM_ival_Stokes, AC_ival, freq_min=freq_min,
                           freq_max=freq_max, dB=True, semilogy=True)

# Repeat this plot focusing on one frequency range
freq_min = 11.5e9  # Hz
freq_max = 13.5e9  # Hz
plotting.plot_gain_spectra(sim_AC, SBS_gain, SBS_gain_PE, SBS_gain_MB, linewidth_Hz,
                           EM_ival_pump, EM_ival_Stokes, AC_ival, freq_min=freq_min,
                           freq_max=freq_max, suffix='_zoom')

print(nbapp.final_report())

4.3.2. Tutorial 3a – Investigating Dispersion and np.save/np.load

This example, contained in tutorials/simo-tut_03_1-dispersion-npload.py calculates the elastic dispersion diagram – the relation between the acoustic wave number \(q\) and frequency \(\Omega\)– for the problem in the previous tutorial. This is done by scanning over the elastic wavenumber q_AC and finding the eigenfrequencies for each value.

As discussed in Formal selection rules for Brillouin scattering in integrated waveguides and structured fibers by C. Wolff, M. J. Steel, and C. G. Poulton DOI:/10.1364/OE.22.032489, the elastic modes of any waveguide may be classified according to their representation of the point group symmetry class corresponding to the waveguide profile. For this problem, the waveguide is rectangular with symmetry group \(C_{2v}\) which has four symmetry classes, which are marked in the dispersion diagram.

This example also takes advantage of the ability to load and save simulation results to save repeated calculation. using the save_simulation and load_simulation methods defined in the mode_calcs module. The previous tutorial saved its electromagnetic results in the file tut02_wguide_data.npz using the Simulation.save_simulation() method, while these tutorial has recovered those results using mode_calc.load_simulation(). This can be a very useful technique when trying to adjust the appearance of plots without having to repeat the whole calculation effort.

Note: from now on, we do not include the code for each tutorial and refer the reader to the relevant files in the <NumBAT>/tutorials directory.

Acoustic dispersion diagram with modes categorised by symmetry as in Table 1 of Wolff et al. Opt. Express. 22, 32489 (2014).

4.3.3. Tutorial 3b – Investigating Dispersion and Multiprocessing

This tutorial, contained in simo-tut_03_2-dispersion-multicore.py continues the study of acoustic dispersion and demonstrates the use of Python multiprocessor calls using the multiprocessing library to increase speed of execution.

In this code as in the previous example, the acoustic modal problem is repeatedly solved at a range of different \(q\) values to build up a set of dispersion curves \(\nu_m(q)\). Due to the large number of avoided and non-avoided crossings, it is usually best to plot dispersion curves like this with dots rather than joined lines. The plot generated below can be improved by increasing the number of \(q\) points sampled through the value of the variable n_qs, limited only by your patience.

The multiprocessing library runs each task as a completely separate process on the computer. Depending on the nature and number of your CPU, this may improve the performance considerably. This can also be easily extended to multiple node systems which will certainly improve performance. A very similar procedure using the threading library allows the different tasks to run as separate threads within the one process. However, due to the existence of the Python Global Interpreter Lock (GIL) which constrains what kinds of operations may run in parallel within Python , multiple threads will typically not improve the performance of NumBAT.

This tutorial also shows an example of saving data, in this case the array of acoustic wavenumbers and frequencies, to a text file using the numpy routine np.savetxt for later analysis.

Acoustic dispersion diagram. The elastic wave number \(q\) is scaled by the phase-matched SBS wavenumber \(2\beta\) where \(\beta\) is the propagation constant of the optical pump mode.

4.3.4. Tutorial 4 – Parameter Scan of Widths

This tutorial, contained in simo-tut_04_scan_widths.py demonstrates the use of a parameter scan of a waveguide property, in this case over the width of the silicon rectangular waveguide, to characterise the behaviour of the Brillouin gain.

The results are displayed in a 3D plot. This may not be the most effective approach for this small data set but gives a sense of what is possible graphically. For a more effective plot, you might like to try the same calculation with around 30 values for the width rather than just 6.

Gain spectra as function of waveguide width.

4.3.5. Tutorial 5 – Convergence Study

This tutorial, contained in simo-tut_05_convergence_study.py demonstrates a scan of numerical parameters for our by now familiar silicon-in-air problem to test the convergence of the calculation results. This is done by scanning the value of the lc_refine parameters. The number of mesh elements (and simulation time) increases with roughly the square of the mesh refinement factor.

For the purpose of convergence estimates, the values calculated at the finest mesh (the rightmost values) are taken as the exact values, notated with the subscript 0, eg. \(\beta_0\). The graphs below show both relative errors and absolute values for each quantity.

Once the convergence properties for a particular problem have been established, it can be useful to do exploratory work more quickly by adopting a somewhat coarser mesh, and then increase the resolution once again towards the end of the project to validate results before reporting them.

Convergence of relative (blue) and absolute (red) optical wavenumbers \(k_{z,i}\). The left axis displays the relative error \((k_{z,i}-k_{z,0})/k_{z,0}\). The right axis shows the absolute values of \(k_{z,i}\).

Convergence of relative (solid, left) and absolute (chain, right) elastic mode frequencies \(\nu_{i}\).

Convergence of photoelastic gain \(G^\text{PE}\). The absolute gain on the right hand side increases down the page because of the convention that NumBAT associates backward SBS with negative gain.

Absolute and relative convergence of moving boundary gain \(G^\text{MB}\).

Absolute and relative convergence of total gain \(G\).

4.3.6. Tutorial 6 – Silica Nanowire

In this tutorial, contained in simo-tut_06_silica_nanowire.py we start to explore the Brillouin gain properties in a range of different structures, in this case a silica nanowire surrounded by vacuum.

The gain-spectra plot below shows the Brillouin gain as a function of Stokes shift. Each resonance peak is marked with the number of the acoustic mode associated with the resonance. This is very helpful in identifying which acoustic mode profiles to examine more closely. In this case, modes 5, 8 and 23 give the most signifcant Brillouin gain. The number of modes labelled in the gain spectrum can be controlled using the parameter mark_mode_thresh in the function plotting.plot_gain_spectra to avoid many labels from modes giving negligible gain. Other parameters allow selecting only one type of gain (PE or MB), changing the frequency range, and plotting on log or dB scales.

It is important to remember that the total gain is not the simple sum of the photoelastic (PE) and moving boundary (MB) gains. Rather it is the coupling terms \(Q_\text{PE}\) and \(Q_\text{MB}\) which are added before squaring to give the total gain. Indeed the two effects may have opposite sign giving net gains smaller than either contribution.

Gain spectrum showing the gain due to the photoelastic effect (PE), the moving boundary effect (PB), and the net gain (Total).

Electromagnetic mode profile of the pump and Stokes field in the \(x\)-polarised fundamental mode of the waveguide.

Mode profiles for acoustic mode 5 which is visible as a MB-dominated peak in the gain spectrum.

Mode profiles for acoustic mode 8 which is visible as a PE-dominated peak in the gain spectrum.

4.3.7. Tutorial 7 – Slot Waveguide

This tutorial, contained in simo-tut_07-slot.py examines backward SBS in a more complex structure: chalcogenide soft glass (\(\text{As}_2\text{S}_3\)) embedded in a silicon slot waveguide on a silica slab. This structure takes advantage of the slot effect which expels the optical field into the lower index medium, enhancing the fraction of the EM field inside the soft chalcogenide glass which guides the acoustic mode and increasing the gain.

Comparing the \(m=2\) and \(m=5\) acoustic mode profiles with the pump EM profile, it is apparent that the field overlap is favourable, where as the \(m=1\) mode gives zero gain due to its anti-symmetry relative to the pump field.

Gain spectrum showing the gain due to the photoelastic effect (PE), the moving boundary effect (PB), and the net gain (Total).

Electromagnetic mode profile of the pump and Stokes field.

Acoustic mode profiles for mode 0.

Acoustic mode profiles for mode 2.

Acoustic mode profiles for mode 1.

Acoustic mode profiles for mode 5.

4.3.8. Tutorial 8 – Slot Waveguide Cover Width Scan

This tutorial, contained in simo-tut_08-slot_coated-scan.py continues the study of the previous slot waveguide, by examining the dependence of the acoustic spectrum on the thickness of a silica capping layer. As before, this parameter scan is accelerated by the use of multi-processing.

It is interesting to look at different mode profiles and try to understand why the eigenfrequency of some modes are more affected by the capping layer. The lowest mode, for instance, is noticeably unaffected.

Acoustic frequencies as function of covering layer thickness.

Modal profiles of lowest acoustic mode.

Modal profiles of second acoustic mode.

Modal profiles of third acoustic mode.

4.3.9. Tutorial 9 – Anisotropic Elastic Materials

This tutorial, contained in simo-tut_09-anisotropy.py improves the treatment of the silicon rectangular waveguide by accounting for the anisotropic elastic properties of silicon (simply by referencing a different material file for silicon).

The data below compares the gain spectrum compared to that found with the simpler isotropic stiffness model used in Tutorial 2. The results are very similar but the isotropic model shows two smaller peaks at high frequency.

Gain spectrum with anisotropic stiffness model of silicon.

Gain spectrum from Tutorial 2 with isotropic stiffness model of silicon.

4.3.10. Tutorial 11 – Two-layered ‘Onion’

This tutorial, contained in simo-tut_11a-onion2.py demonstrates use of a two-layer onion structure for a backward intra-modal SBS configuration. Note that with the inclusion of the background layer, the two-layer onion effectively creates a three-layer geometry with core, cladding, and background surroundings. This is the ideal structure for investigating the cladding modes of an optical fibre. It can be seen by looking through the optical mode profiles in tut_11a-fields/EM*.png that this particular structure supports five cladding modes in addition to the three guided modes (the \(\mbox{TM}_0\) mode is very close to cutoff).

Next, the gain spectrum and the mode profiles of the main peaks indicate as expected, that the gain is optimal for modes that are predominantly longitudinal in character.

The accompanying tutorial simo-tut_11b-onion3.py introduces one additional layer and would be suitable for studying the influence of the outer polymer coating of an optical fibre or depressed cladding W fibre.

Gain spectrum for the two-layer structure in tut_11a.

Mode profile for fundamental optical mode.

Mode profile for acoustic mode 17.

Mode profile for acoustic mode 32.

Mode profile for acoustic mode 39.

Gain spectrum for the three-layer structure in tut_11b.

Mode profile for acoustic mode 44.

Mode profile for acoustic mode 66.

Mode profile for acoustic mode 66.

4.3.11. Tutorial 12 – Valdiating the calculation of the EM dispersion of a two-layer fibre

How can we be confident that NumBAT’s calculations are actually correct? This tutorial and the next one look at rigorously validating some of the modal calculations produced by NumBAT.

This tutorial, contained in simo-tut_12.py, compares analytic and numerical calculations for the dispersion relation of the electromagnetic modes of a cylindrical waveguide. This can be done in both a low-index contrast (SMF-28 fibre) and high-index contrast (silicon rod in silica) context. We calculate the effective index \(\bar{n}\) and normalised waveguide parameter \(b=(\bar{n}^2-n_\text{cl}^2)/(n_\text{co}^2-n_\text{cl}^2)\) as a function of the normalised freqency \(V=k a \sqrt{n_\text{co}^2-n_\text{cl}^2}\) for radius \(a\) and wavenumber \(k=2\pi/\lambda\). As in several previous examples, this is accomplished by a scan over the wavenumber \(k\).

The numerical results (marked with crosses) are compared to the modes found from the roots of the rigorous analytic dispersion relation (solid lines). We also show the predictions for the group index \(n_g=\bar{n} + V \frac{d\bar{n}}{dV}\). The only noticeable departures are right at the low \(V\)-number regime where the fields become very extended into the cladding and interact significantly with the boundary. The results could be improved in this regime by choosing a larger domain at the expense of a longer calculation.

As this example involves the same calculation at many values of the wavenumber \(k\), we again use parallel processing techniques. However, in this case we demonstrate the use of threads (multiple simultaneous strands of execution within the same process) rather than a pool of separate processes. Threads are light-weight and can be started more efficiently than separate processes. However, as all threads share the same memory space, some care is needed to prevent two threads reading or writing to the same data structure simultaneously. This is dealt with using the helper functions and class in the numbattools.py module.

Optical effective index as a function of normalised frequency for SMF-28 fibre.

Optical normalised waveguide parameter as a function of normalised frequency for SMF-28 fibre.

Optical group index \(n_g=\bar{n} + V \frac{d\bar{n}}{dV}\) as a function of normalised frequency for SMF-28 fibre.

Optical effective index as a function of normalised frequency for silicon rod in silica.

Optical normalised waveguide parameter as a function of normalised frequency for silicon rod in silica.

Optical group index \(n_g=\bar{n} + V \frac{d\bar{n}}{dV}\) as a function of normalised frequency for silicon rod in silica.

4.3.12. Tutorial 13 – Valdiating the calculation of the dispersion of an elastic rod in vacuum

The tutorial simo-tut_13.py performs the same kind of calculation as in the previous tutorial for the acoustic problem. In this case there is no simple analytic solution possible for the two-layer cylinder. Instead we create a structure of a single elastic rod surrounded by vacuum. NumBAT removes the vacuum region and imposes a free boundary condition at the boundary of the rod. The modes found are then compared to the analytic solution of a free homonegenous cylindrical rod in vacuum.

We find excellent agreement between the analytical (coloured lines) and numerical (crosses) results. Observe the existence of two classes of modes with azimuthal index \(p=0\), corresponding to the pure torsional modes, which for the lowest band propagate at the bulk shear velocity, and the so-called Pochammer hybrid modes, which are predominantly longitudinal, but must necessarily involve some shear motion to satisfy mass conservation.

It is instructive to examine the mode profiles in tut_13-fields and track the different field profiles and degeneracies found for each value of \(p\). By basic group theory arguments, we know that every mode with \(p\ne 0\) must come as a degenerate pair and this is satisfied to around 5 significant figures in the calculated results. It is interesting to repeat the calculation with a silicon (cubic symmetry) rod rather than chalcogenide (isotropic). In that case, the lower symmetry of silicon causes splitting of a number of modes, so that a larger number of modes are found to be singly degenerate, though invariably with a partner state at a nearby frequency.

Elastic frequency as a function of normalised wavenumber for a chalcogenide rod in vacuum.

Elastic “effective index” defined as the ratio of the bulk shear velocity to the phase velocity \(n_\text{eff}=V_s/V\), for a chalcogenide rod in vacuum.

4.3.13. Tutorial 14 – Multilayered ‘Onion’

This tutorial, contained in simo-tut_14-multilayer-fibre.py shows how one can create a circular waveguide with many concentric layers of different thickness. In this case, the layers are chosen to create a concentric Bragg grating of alternating high and low index layers. As shown in C. M. de Sterke, I. M. Bassett and A. G. Street, “Differential losses in Bragg fibers”, J. Appl. Phys. 76, 680 (1994), the fundamental mode of such a fibre is the fully azimuthally symmetric \(\mathrm{TE}_0\) mode rather than the usual \(\mathrm{HE}_{11}\) quasi-linearly polarised mode that is found in standard two-layer fibres.

Fundamental electromagnetic mode profile for the concentric Bragg fibre.

4.3.14. Tutorial 15 – Coupled waveguides

This tutorial, contained in tutorials/simo-tut_14-coupled-wg.py demonstrates the supermode behaviour of both electromagnetic and elastic modes in a pair of closely adjacent waveguides.

4.3.15. Tutorial 16 – Using NumBAT in a Jupyter notebook

For those who like to work in an interactive fashion, NumBAT works perfectly well inside a Jupyter notebook. This is demonstrated in the file jup_16_smf28.ipynb using the standard SMF-28 fibre problem as an example.

On a Linux system, you can open this at the command line with:

$ jupyter-notebook ``jup_16_smf28.ipynb``

or else load it directly in an already open Jupyter environment.

The notebook demonstrates how to run standard NumBAT routines step by step. The output is still written to disk, so the notebook includes some simple techniques for efficiently displaying mode profiles and spectra inside the notebook.

NumBAT calculation of SMF-28 modes within a Jupyter notebook.