Tutorial

Installation

TransferMatrix.jl is a part of Julia's general registry and the source code can be found at https://github.com/garrekstemo/TransferMatrix.jl. From the Julia REPL, enter the package manager mode mode by typing ]. Then just enter the following to install the package:

pkg> add TransferMatrix

A simple calculation

To get you up and running, let's build a simple two-layer structure of air and glass and calculate the reflectance and transmittance to visualize the Brewster angle for p-polarized light. We fix the wavelength of incident light and vary the angle of incidence when setting up our Structure. (It is just as simple to fix the incidence angle and calculate the transfer matrix as a function of the field wavelength.)

We start by making a Layer type of air and a Layer of glass. We'll do this for a wavelength of 1 μm. Since there are only two layers and the transfer matrix method treats the first and last layers as semi-infinite, there is no need to provide a thickness for our glass and air layers. From the examples below, you can see that there are fields for

the material name
the layer thickness
a list of wavelengths
the real part of the refractive index (corresponding to the wavelength)
the imaginary part of the refractive index

Details about different ways to make a layer are further on in the tutorial.

using TransferMatrix

air = Layer("Air", 0.0, [1.0e-6], [1.0], [0.0])
glass = Layer("Glass", 0.0, [1.0e-6], [1.5], [0.0])

Layer("Glass", 0.0, [1.0e-6], [1.5], [0.0])

Assembling layers into a structure

Now that we have our glass and air layers, we need to assemble them into a structure and provide the angles of the field with respect to the surface of the structure. We do this with the Structure type.

θs = range(0., 85., length = 500)
s = Structure([air, glass], [1.0e-6], collect(θs) .* π/180)

Structure(Layer[Layer("Air", 0.0, [1.0e-6], [1.0], [0.0]), Layer("Glass", 0.0, [1.0e-6], [1.5], [0.0])], [1.0e-6], [0.0, 0.0029730057398701004, 0.005946011479740201, 0.0089190172196103, 0.011892022959480402, 0.0148650286993505, 0.0178380344392206, 0.020811040179090703, 0.023784045918960803, 0.026757051658830907  …  1.4567728125363493, 1.459745818276219, 1.4627188240160895, 1.4656918297559594, 1.4686648354958296, 1.4716378412356999, 1.4746108469755699, 1.47758385271544, 1.48055685845531, 1.4835298641951802])

The first argument is just a list of layers. The second argument is a list of desired wavelengths. Often the refractive index data we have for two materials are not given for exactly the same wavelengths. TransferMatrix.jl uses an interpolation function to normalize the wavelengths and complex refractive indices for all layers from this user-provided list of wavelengths. (Be careful if the range you provide goes beyond the range of the data that you have!) Now we can evaluate the angle-resolved spectrum using the function angle_resolved().

res = angle_resolved(s)

TransferMatrix.AngleResolvedResult([0.04000000000000003; 0.039999528598952526; … ; 0.4813707135624852; 0.4932538118188537;;], [0.04000000000000003; 0.04000047140359379; … ; 0.7246514855242108; 0.7323454787108021;;], [0.9600000000000001; 0.9600004714010474; … ; 0.5186292864375149; 0.5067461881811459;;], [0.9600000000000001; 0.9599995285964057; … ; 0.2753485144757894; 0.267654521289198;;], Matrix{ComplexF64}[], ComplexF64[0.0 + 0.0im; 0.00297300136025645 + 0.0im; … ; 0.9959311813887618 + 0.0im; 0.9961946980917455 + 0.0im;;])

Let's also plot the result using the Makie.jl data visualization package.

using CairoMakie

brewster = atan(1.5) * 180 / π

f = Figure()
ax = Axis(f[1, 1], xlabel = "Incidence Angle (°)", ylabel = "Reflectance / Transmittance")

lines!(θs, res.Tss[:, 1], color = :firebrick4, label = "Ts")
lines!(θs, res.Tpp[:, 1], color = :orangered3, label = "Tp")
lines!(θs, res.Rss[:, 1], color = :dodgerblue4, label = "Rs")
lines!(θs, res.Rpp[:, 1], color = :dodgerblue1, label = "Rp")
vlines!(brewster, color = :dodgerblue1, linestyle = :dash)
text!("Brewster angle\n(Rp = 0)", position = (35, 0.6))

axislegend(ax)
f

We can see that the result of the angle-resolved calculation has four solutions: the s-wave and p-wave for both the reflected and transmitted waves. And we see that the Brewster angle is $\arctan\left( n_\text{glass} /n_\text{air} \right) \approx 56^{\circ}$, as expected. Simultaneous calculation of s- and p-polarized incident waves is a feature of the general 4x4 transfer matrix method being used. The angle_resolved function will also loop through all wavelengths so that you can plot a color plot of wavelength and angle versus transmittance (or reflectance).

Defining a Layer

The Layer type is immutable. Once you make one, you can't change any of its characteristics later. There are two ways to define a Layer. The first is by directly filling in the argument fields:

Layer(material, thickness, λs, ns, κs)

You may choose any units you like for thickness and wavelength, but they must be the same (e.g. both are in nanometers). The material is just a String. Wavelength, and the refractive index arguments must be Arrays, even if they just contain with a single item. A very manual way to make a Layer might look like this:

glass = Layer("Glass", 20.0e-6, [1.0e-6, 1.1e-6, 1.3e-6], [0.0, 0.0, 0.0])

This way works well for simple Layers or when you just need a single frequency for an angle-resolved calculation, but this is a lot more work if the refractive index is pulled from a database or a file. The website refractiveindex.info contains a large database of refractive indices from peer-reviewed papers. TransferMatrix.jl makes it easy to use a CSV file downloaded from this website with the read_refractive() function, which returns a Layer type. You can choose to include (or not) the real/imaginary part of the complex refractive index. Note that this refractiveindex.info stores wavelength information in units of micrometers.

glass = read_refractive("/path/to/file.csv", "Glass", 20e-6, div=1e6)

The div keyword argument at the end of read_refractive represents the desired unit conversion. In this case, if we want everything to be in meters we must divide each wavelength in the raw data by $10^6$. Alternatively, we can change the thickness units to micrometers. Just be careful that all Layers in your multi-layered structure have consistent units.

The Layer properties can be independently queried. To get the material name or thickness, for example, you can just type:

julia> glass.material
"Glass"
julia> glass.thickness
2.0e-5

The same can be done to get wavelengths (glass.λ) (get λ by typing \lambda), real refractive index (glass.n), and complex refractive index (glass.κ) (get κ by typing \kappa). (A very useful feature of Julia is built-in support for unicode characters).

A simple multi-layered structure

Now that we can make Layers, we can put them together to calculate the transfer matrix for a multi-layered structure and plot the reflectance and transmittance spectra. A very important structure in optics is the Fabry-Pérot etalon, which can be made with two mirrors that face each other with a gap between them. Let's make this a little more challenging (and experimentally realistic) by mounting the mirrors onto windows that are transparent in our chosen wavelength range. The gap will simply be 20 μm of air. We will scan in the mid infrared between 4 and 6 μm and use data generated via the Lorentz-Drude model for each 10 nm-thick gold layer and experimental data from Malitson, 1963 for the windows. (Some example refractive index data are in the refractive_index_data folder in the simulation src directory. The ../ syntax denotes a relative file path one level up from the current folder.)

λs = range(4.0, 6.0, length = 1000) .* 1e-6
aufile = "../../../../refractive_index_data/Au_nk_0.248-6.20um_Lorentz-Drude_Rakic1998.csv"
caf2file = "../../../../refractive_index_data/CaF2_n_0.23-9.7um_Malitson.csv"

air = Layer("Air", 20e-6, collect(λs), fill(1.0, length(λs)), zeros(length(λs)))
au = read_refractive(aufile, "Au", 10e-9, div=1e6)
caf2 = read_refractive(caf2file, "CaF2", 5.0e-6, div=1e6)

s = Structure([caf2, au, air, au, caf2], collect(λs), [0.0])
Tp, Ts, Rp, Rs = calculate_tr(s)

f, ax, l = lines(λs .* 1e6, Ts)
ax.xlabel = "Wavelength"
ax.ylabel = "Transmittance"
f

It is easy to combine manually-generated data and experimental data to calculate the global transfer matrix.

Periodic structures with a YAML config file

With the tools described above, it is pretty easy to make any kind of multi-layered structure that you want. A very large structure, maybe one with repeating layers, is also not difficult to build. However, it would be useful to be able to load a structure from a set of parameters so that it can be saved or even shared with others for better reproducibility. To do this, a Structure can be loaded directly from a YAML file (using YAML.jl).

Let's demonstrate this with a quarter-wave stack, which is a periodic structure with two alternative layers where the thickness is one fourth of the wavelength within the medium, or the relation

\[d_i = \frac{\lambda}{4n_i},\]

where $d_i$ is the thickness for layer $i$, $\lambda$ is the electric field wavelength in vacuum, and $n_i$ is the index of refraction for layer $i$.

Here is the YAML file that gives us this:

# All units in micrometers

min_wavelength: 0.5
max_wavelength: 2.0
n_wavelengths: 1000

# Theta in radians
theta_i: 0.0
theta_f: 0.0
n_theta: 1

layers:

    layer1:
        material: Air
        thickness: 0.5
        wavelength: 1.0
        refractive_index: 1.0
        extinction_coeff: 0.0

    layer2:
        periods: 3
        layers:
          layer1:
            material: "ZnS n=2.32"
            thickness: 0.1077
            refractive_filename: ZnS_n_0.4-14um_Amotchkina2020.csv
          layer2:
            material: "MgF2 n=1.36"
            thickness: 0.1838
            refractive_filename: MgF2_nk_0.03-2.0um_Rodriguez2017.csv

    layer3: 
        material: "ZnS n=2.32"
        thickness: 0.1838
        refractive_filename: ZnS_n_0.4-14um_Amotchkina2020.csv

At the top, we set the minimum and maximum wavelengths to be calculated, with the total number of points set to 1000. We won't do an angle-resolved calculation, so initial and final $\theta$ are zero (with n_theta just one). Next the multi-layered structure is defined under layers:. The first layer is just air, so we don't need a file for this and can set the refractive index for a single wavelength. These values and wavelength are expanded when the YAML file is loaded. The final layer (layer3) is zinc sulfide (ZnS), which is loaded from a CSV file downloaded from refractiveindex.info. The middle section (layer2) is where we define the periodic structure. The first item determines how many periods we want. In this case, periods: 3. Then we set the layers that alternate under another layers section. You can see that the first layer layer1 is zinc sulfide and the second layer layer2 is magnesium fluoride (MgF₂). The number of times that these two layers are repeated can easily be changed by editing the periods field. You can even have multiple repeating structures within the global structure for arbitrarily complicated structures.

To load the YAML config file into a Structure, we use the load_from_yaml() function. You can find this example in the default_config folder in the quarter-wave.yaml file.

s = load_from_yaml("../../../default_config/quarter-wave.yaml", 1e-6)
Tp, Ts, Rp, Rs = calculate_tr(s)

([0.9130563098195978, 0.9029002689472558, 0.8926804707814351, 0.8826271417027922, 0.8729449221143134, 0.863814941698935, 0.8553728415067637, 0.8477399728437076, 0.8410130207024722, 0.8352772978673564  …  0.8580061039350702, 0.8588749011786784, 0.859743674915508, 0.8606123664190505, 0.8614809171863055, 0.8623492470576426, 0.863217315621725, 0.8640850689508629, 0.8649524494734356, 0.8658193998566369], [0.9130563098195978, 0.9029002689472558, 0.8926804707814351, 0.8826271417027922, 0.8729449221143134, 0.863814941698935, 0.8553728415067637, 0.8477399728437076, 0.8410130207024722, 0.8352772978673564  …  0.8580061039350702, 0.8588749011786784, 0.859743674915508, 0.8606123664190505, 0.8614809171863055, 0.8623492470576426, 0.863217315621725, 0.8640850689508629, 0.8649524494734356, 0.8658193998566369], [0.07738850881115539, 0.08767249740770969, 0.09802099094735575, 0.10820161347294284, 0.1180088042284185, 0.12726665774657298, 0.13582461972538507, 0.1435627549446024, 0.1503868114592622, 0.15621076026724637  …  0.14157676668653812, 0.14070839627869383, 0.13984004986029905, 0.13897178617577158, 0.13810366374575572, 0.13723576150392724, 0.13636812090021944, 0.1355007961012557, 0.13463384469529158, 0.13376732403151376], [0.07738850881115539, 0.08767249740770969, 0.09802099094735575, 0.10820161347294284, 0.1180088042284185, 0.12726665774657298, 0.13582461972538507, 0.1435627549446024, 0.1503868114592622, 0.15621076026724637  …  0.14157676668653812, 0.14070839627869383, 0.13984004986029905, 0.13897178617577155, 0.13810366374575572, 0.13723576150392724, 0.13636812090021944, 0.1355007961012557, 0.13463384469529158, 0.13376732403151376])

This particular example is taken from Pochi Yeh's Optical Waves in Layered Media on page 110. You can find the reflectivity of this structure in Table 5.1. The values we have calculated are slightly different since the data used here is slightly dispersive with wavelength, but the example in the book takes flat refractive index values. You can try plotting this structure for an increasing number of periods and observe how the reflectance near $\lambda$ = 1 μm increases.

using CairoMakie

f = Figure()
ax = Axis(f[1, 1], title = "ZnS / MgF quarter-wave with 3 layers", xlabel = "Wavelength (nm)", ylabel = "Transmittance / Reflectance")

lines!(s.λ .* 1e9, Tp, label = "T")
lines!(s.λ .* 1e9, Rp, label = "R")
axislegend(ax, position = :rc)

f

Electric field

The electric field can be calculated as a function of position within the layered structure using the electric_field() function, which takes the Structure and desired wavelength λ, as well as optional argument angle of incidence θ and optional keyword argument for the number of data points numpoints. We can plot the field profile for the distributed bragg reflector (DBR) we constructed in the previous section. Let's do this for λ = 1 μm.

λ_field = 1e-6
field = electric_field(s, λ_field)

f = Figure()
ax = Axis(f[1, 1], title = "Electric Field Profile at λ = $(Int(λ_field * 1e9)) nm", xlabel = "z position (nm)", ylabel = "Field intensity (a.u.)")

lines!(field.z .* 1e9, real(field.p[1, :]).^2)

vlines!(field.boundaries[1], color = :gray30, linestyle = :dash)
vlines!(field.boundaries[end] * 1e9, color = :gray30, linestyle = :dash)

f

The electric field result contains the position z within the structure, the (x, y, z) components (corresponding to the first, second, and third components of the Array) of the p-polarized and s-polarized light, and the positions of the layer interfaces. So for example, to get the x-component of the p-polarized field along all of z, we would call field.p[1, :], as we have done above.

The layer boundaries is useful for plotting (as shown in the figure above) and checking that the in-plane components are continuous throughout the structure, as required by Maxwell's interface conditions.

Minor implementation details

Behind the scenes of the electric field calculation, a new Structure is actually being created and initialized for the user-supplied wavelength. The transfer matrix is again calculated for the provided wavelength and light incidence angle. Some solution of this sort is necessary because the exact wavelength that you may wish to calculate for may not be part of the original Structure.