1. Introduction to NumBAT

1.1. Introduction

NumBAT, the Numerical Brillouin Analysis Tool, integrates electromagnetic and acoustic mode solvers to calculate the interactions of optical and acoustic waves in waveguides.

1.2. Goals

NumBAT is designed primarily to calculate the optical gain response from stimulated Brillouin scattering (SBS) in integrated waveguides. It uses finite element algorithms to solve the electromagnetic and acoustic modes of a wide range of 2D waveguide structures. It can account for photoelastic/electrostriction and moving boundary/radiation pressure effects, as well as arbitrary acoustic anisotropy.

NumBAT also supports user-defined material properties and we hope its creation will drive a community-driven set of standard properties and geometries which will allow all groups to test and validate each other’s work.

A full description of the NumBAT physics and numerical algorithms is available in the article B.C.P Sturmberg at al., “Finite element analysis of stimulated Brillouin scattering in integrated photonic waveguides”, J. Lightwave Technol. 37, 3791-3804 (2019), available at https://dx.doi.org/10.1109/JLT.2019.2920844.

NumBAT is open-source software and the authors welcome additions to the code. Details for how to contribute are available in Contributing to NumBAT.

1.3. Citing NumBAT

If you use NumBAT in published work, we would appreciate a citation to B.C.P Sturmberg at al., “Finite element analysis of stimulated Brillouin scattering in integrated photonic waveguides”, J. Lightwave Technol. 37, 3791-3804 (2019), available at https://dx.doi.org/10.1109/JLT.2019.2920844 and https://arxiv.org/abs/1811.10219, and a link to the github page at https://github.com/michaeljsteel/NumBAT.

1.4. Development team

NumBAT was developed by Bjorn Sturmberg, Kokou Dossou, Blair Morrison, Chris Poulton and Michael Steel in a collaboration between Macquarie University, the University of Technology Sydney, and the University of Sydney.

We thank Christian Wolff, Mike Smith and Mikolaj Schmidt for contributions.

1.5. Contributing to NumBAT

NumBAT is open source software licensed under the GPL with all source and documentation available at github.com. We welcome additions to NumBAT code, documentation and the materials library. Interested users should fork the standard release from github and make a pull request when ready. For major changes, we strongly suggest contacting the NumBAT team before starting work at michael.steel@mq.edu.au.

1.6. Support

Development of NumBAT has been supported in part by the Australian Research Council under Discovery Projects DP130100832, DP160101691, DP200101893 and DP220100488.

1.7. Release notes

1.7.1. Version 2.0

A number of API changes have been made in NumBAT 2.0 to tidy up the interface and make plotting and analysis simpler and more powerful. You will need to make some changes to existing files to run in NumBAT 2.0. Your best guide to new capabilities and API changes is to look through the code in the tutorial examples.

Some key changes you will need to make are as follows:

• On Linux, the fortran Makefile is now designed to work with a virtual environment python to avoid dependencies on your system python.

• There is a new core NumBAT module numbat that should be imported before any other NumBAT modules.

• It should no longer be necessary to import the object or Numbat (note different case) modules.

• The first call to any NumBAT code should be to create a NumBAT application object by calling nbapp = numbat.NumBATApp().

• The default output prefix can now be set as an argument to numbat.NumBATApp(). All output can be directed to a sub-folder of the starting directory with a second argument: nbapp = numbat.NumBATApp('tmp', 'tmpdir').

• The waveguide class Struct has been renamed to Structure.

• A waveguide is now constructed using nbapp.make_waveguide rather than object.Structure.

• The interface for creating materials has changed. You now call the materials. make_material( name ) function. For example material_a = materials.make_material('Vacuum')

• To access an existing material in an existing Struture object (say, in a variable called wguide) use wguide.get_material( label ) For example, mat_a = wguide.get_material('b') where the allowed labels are bkg and the letters a to r.

• The member name for refractive index in a Material object has changed from n to refindex_n.

• The member name for density in a Material object has changed from n to rho.

• Due to a change in parameters, the function plotting.gain_spectra is deprecated and replaced by plotting.plot_gain_spectra with the following changes:

  • The frequency arguments freq_min and freq_max should now be passed in units of Hz, not GHz.

  • The argument k_AC has been removed.

• In all functions the parameter prefix_str has been renamed to prefix for brevity. Using the default output settings in NumBATApp(), these should be rarely needed.

• All waveguides are now specified as individual plugin classes in the files backend/msh/user_waveguides.json and backend/msh/user_meshes.py. These files provide useful examples of how to design and load new waveguide templates. See the following chapter for more details.

1.8. What does NumBAT actually calculate?

NumBAT performs three main types of calculations given a particular waveguide design:

• solve the electromagnetic modal problem using the finite element method (FEM).

• solve the elastic modal problem using FEM.

• calculate Brillouin gain coefficients and linewidths for a given triplet of two optical and one elastic mode, and use this to generate gain spectra.

Here we specify the precise mathematical problems been solved. For further details, see the NumBAT paper in the Journal of Lightwave Technolgoy at at https://dx.doi.org/10.1109/JLT.2019.2920844.

1.8.1. Electromagnetic modal problem

The electromagnetic wave problem is defined by the vector wave equation

\[- \nabla \times (\nabla \times {\vec E}) + \omega^2 \epsilon_0 \, \epsilon_r(x,y) \vec E =0,\]

where the electric field has the form for modal propagation along \(z\):

\[\vec E(x,y,z,t) = a(z) \vec e(x,y) e^{i (kz-\omega t) } + a^*(z) \vec e(x,y) e^{-i (kz-\omega t) } ,\]

and the magnetic field is given by

\[\vec B = \frac{1}{i \omega} \nabla \times \vec E.\]

1.8.2. Elastic modal problem

The elastic modal problem is defined by the wave equation

\[\nabla \cdot \bar{T} + \omega^2 \rho(x,y) \vec U = 0,\]

where \(\vec u\) is the elastic displacement and \(\bar{T}=\mathbf{c}(x,y) \bar{S}\) is the stress tensor, defined in terms of the stiffness \(\mathbf{c}\) and the strain tensor \(\bar{S}=S_{ij} = \frac{1}{2}(\frac{\partial U_i}{\partial r_j} + \frac{\partial U_j}{\partial r_i})\).

The displacement has the modal propagation form

\[\vec U = b(z) \, \vec u(x,y) e^{i (qz-\Omega t) } + b^*(z) \, \vec u(x,y) e^{-i (qz-\Omega t) } ,\]

For details on how these problems are framed as finite element problems, we refer to https://dx.doi.org/10.1109/JLT.2019.2920844.

1.8.3. Modal properties

For propagation in a given mode \(\vec e_n\) or \(\vec U_n\), the optical (\(o\)) and elastic (\(a\)) energy fluxes in Watts and linear energy densities in (J/m) are given by the following expressions

\[\begin{split}{\cal P}_n^{(o)} & = 2 \mathrm{Re} \int \mathrm{d}^2 r \, \hat{z} \cdot (\vec e_n^*(x,y) \times \vec h_n(x,y)), \\ {\cal E}_n^{(o)} & = 2 \epsilon_0 \int \mathrm{d}^2 r \, \epsilon_r(x,y) |\vec e_n(x,y)|^2 \\ {\cal P}_n^{(a)} & = \mathrm{Re} \int \mathrm{d}^2 r \, (-2i\Omega) \sum_{jkl} c_{zjkl}(x,y) u^*_{mj}(x,y) \partial_k u_{ml}(x,y) \\ {\cal E}_n^{(a)} & = 2 \Omega^2 \int_A \mathrm{d}^2 r \, \rho(x,y) |\vec u_n(x,y)|^2.\end{split}\]

For fields with slowly-varying optical and elastic amplitudes \(a_m(z)\) and \(b_m(z)\), the total carried powers are

\[\begin{split}P^{(o)}(z) & = \sum_m |a_n(z)|^2 {\cal P}_n^{(o)} \\ P^{(a)}(z) & = \sum_m |b_n(z)|^2 {\cal P}_n^{(a)}.\end{split}\]

Note that in this convention, the amplitude functions \(a_m(z)\) and \(b_m(z)\) are dimensionless and the units of the fields live in the modal functions \(\vec e_n, \vec h_n, \vec U_n\).

1.8.4. SBS gain calculation modal problem

The photoelastic and moving boundary couplings in J/m are given by

\[\begin{split}Q^{\mathrm{(PE)}} & = - \epsilon \int_A \mathrm{d}^2 r \, \sum_{ijkl} \epsilon_r^2 \, e_i^{(s)*} \, e_j^{(p)} \, p_{ijkl} \, \partial_k u_l^* \\ Q^{\mathrm{(MB))}} & = \int_{\cal C} \mathrm{d} {\vec r} \, (\vec u^* \cdot \hat{n}) \times \\ & ~~~~ \left [ (\epsilon_a - \epsilon_b) \epsilon_0 (\hat{n} \times \vec u^{(s)})^* \cdot (\hat{n} \times \vec e^{(p)}) - (\epsilon_a^{-1} - \epsilon_b^{-1}) \epsilon_0^{-1} (\hat{n} \cdot \vec d^{(s)})^* \cdot (\hat{n} \cdot \vec d^{(p)}) \right]\end{split}\]

Then the peak SBS gain \(\Gamma\) is given by

\[\Gamma = \frac{2\omega \Omega}{\alpha_t} \frac{|Q_\mathrm{tot}|^2}{P^{(s)}P^{(p)}{\cal E}^{(a)}},\]

where the total SBS coupling is \(Q_\mathrm{tot} = Q^{(\mathrm{PE})} + Q^{(\mathrm{MB})}\).

Here \(\alpha_t\) is the temporal elastic loss coefficent in \(\mathrm{s}^{-1}\). It is related to the spatial attenuation coefficient by \(\alpha_s = \alpha_t /v_{\mathrm{p}}^{(\mathrm{a})}\) with \(v_{\mathrm{p}}^{(\mathrm{a})}\) being the elastic phase velocity.

In a backward SBS problem, where there is genuine gain in Stokes optical field propagating in the negative \(z\) direction, its power evolves as

\[P^{\mathrm{(s)}}(z) = P_\mathrm{in}^{\mathrm{(s)}} e^{-\Gamma z}.\]

1.8.5. SBS equations of motion

With the above conventions, the dynamical equations for the slowly-varying amplitudes are

\[\begin{split}\frac{1}{v_g^{p}} \frac{\partial }{\partial t} a_p + \frac{\partial }{\partial z} a_p & = -i \frac{\omega_p Q_{\mathrm{tot}}}{{\cal P}_o} a_s b \\ \frac{1}{v_g^{s}} \frac{\partial }{\partial t} a_s - \frac{\partial }{\partial z} a_s & = i \frac{\omega_s Q_{\mathrm{tot}}^*}{{\cal P}_o} a_p b^* \\ \frac{1}{v_g^{a}} \frac{\partial }{\partial t} b + \frac{\partial }{\partial z} b + \frac{\alpha_s}{2} b & = -i \frac{\Omega Q_{\mathrm{tot} }^*}{{\cal P}_a} a_p a_s^*\end{split}\]

Here we’ve chosen the group velocities to be positive and included the propagation direction explicitly.

1.8.6. Connect these to output quantities from code

1.8.7. Equivalent forms of equations

TODO: show forms without the normalisation energies and with cubic style effective area.

Compare to some fiber literature and the hydrodynamic reprn.