# 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`

and are ordered by decreasing effective index and increasing frequency respectively.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\).The integer \(m\) therefore has no particular correspondence to the conventional two index mode indices for fibre or rectangular waveguides.

`inc_a_x, inc_a_y, inc_b_x, inc_b_y, slab_a_x, slab_a_y,`

… etcThese 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:

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:

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).

#### 4.3.1.3. Miscellaneous comments¶

Here are some further elements to note about this example:

When using the

`fast=`

mode, the output data and fields directory begin with`ftut_02`

rather than`tut_02`

.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.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.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

magneticfield \(\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.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.By default, plots are exported as

`png`

format. Pass the option`pdf_png=pdf`

to plot functions to generate aPlots 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.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.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.

### 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.

### 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.

### 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.

### 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.

### 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.

### 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.

### 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.

### 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.

### 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.

### 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.

### 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.

### 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.