Beam propagation

Finesse provides a tool-set of beam tracing and mode matching functions which can be used outside of a simulation context. These functions allow for arbitrary beam propagation over any optical path of a given configuration - supporting both numeric and symbolic calculations.

The primary function recommended for most beam propagation analyses is propagate_beam(), also callable from a model via Model.propagate_beam().

The ‘propagate_beam’ function

As noted in the API documentation linked above, the propagate_beam() function can take any optical path of a defined model and propagate an arbitrary beam from the start node to the end node of the path.

Here we will use a proposed, preliminary design [24] of the Einstein Telescope Low-Frequency (ET-LF) signal recycling cavity (SRC), complete with arm telescopes, to highlight how this function can be used. The model is defined below via Finesse kat-file syntax (see KatScript).

import finesse
finesse.configure(plotting=True)

model = finesse.Model()
model.parse("""
### L0 -> BS -> YARM of ET-LF
# input
l L0 P=1
s l0 L0.p1 BS.p1 L=10

# Main beam splitter
bs BS T=0.5 L=37.5u alpha=60
s BSsub1 BS.p3 BSAR1.p1 L=0.07478 nr=nsilica
s BSsub2 BS.p4 BSAR2.p1 L=0.07478 nr=nsilica
bs BSAR1 R=50u L=0 alpha=-36.6847
bs BSAR2 R=50u L=0 alpha=36.6847

# Y arm telescope
s lBS_ZM1 BS.p2 ZM1.p1 L=70
bs ZM1 T=250u L=37.5u Rc=-50
s lZM1_ZM2 ZM1.p2 ZM2.p1 L=50
bs ZM2 T=0 L=37.5u Rc=-82.5
s lZM2_ITMlens ZM2.p2 ITM_lens.p1 L=52.5

lens ITM_lens 75
s lITM_th2 ITM_lens.p2 ITMAR.p1 L=0

# Y arm input mirror
m ITMAR R=0 L=20u
s ITMsub ITMAR.p2 ITM.p1 L=0.2 nr=nsilicon
m ITM T=7000u L=37.5u Rc=-5580

# Y arm length
s l_arm ITM.p2 ETM.p1 L=10k

# Y arm end mirror
m ETM T=6u L=37.5u Rc=5580
s ETMsub ETM.p2 ETMAR.p1 L=0.2 nr=nsilicon
m ETMAR R=0 L=500u

# SRM
s lBS_SRM BSAR2.p3 SRM.p1 L=10

m SRM T=0.2 L=0 Rc=-9410
s SRMsub SRM.p2 SRMAR.p1 L=0.0749 nr=nsilicon
m SRMAR R=0 L=50n

# cavities
cav cavARM ITM.p2
cav cavSRC SRM.p1 ITM.p1.i

var nsilica 1.44963098985906
var nsilicon 3.42009

lambda(1550n)
""")

We can call propagate_beam() using any two optical nodes as the targets, in this case we propagate the beam from the ITM to the SRM.

ps = model.propagate_beam(model.ITM.p1.o, model.SRM.p1.i)

Note that we can use Port objects as the end points too, in which case the from_node argument is deduced to be the output optical node of the port and to_node will be the input optical node of the other port. In our example, the line below is equivalent to the above.

ps = model.propagate_beam(model.ITM.p1, model.SRM.p1)

By specifying the two end nodes / ports, the optical path between these points will be determined in the function. This is generally a fast operation, however it is likely the slowest of all operations carried out by the propagate_beam() function as a whole. Therefore, in the rare case where this function is being called in a large loop, it makes sense to pre-compute the optical path and pass this instead. An example of this is shown below, again equivalent to the above.

ITM_TO_SRM = model.path(model.ITM.p1, model.SRM.p1)
ps = model.propagate_beam(path=ITM_TO_SRM)

One thing you may note from this type of call to propagate_beam() is that no input beam parameter (q_in) argument has been specified. This means that this argument will be automatically deduced from the model via an internal Model.beam_trace() call - where the beam parameter at the input node of the path is then accessed and used as the q_in argument. A custom beam parameter can be used by specifying a q_in argument explicity (this can be a complex number or a BeamParam instance).

Another aspect to take into account is that propagate_beam() operates on a single plane (i.e. ‘x’ for the tangential plane, ‘y’ for the sagittal plane). By default, the tangential plane is used. This can be changed to sagittal by specifying direction='y' as an argument.

The return value of propagate_beam() is a PropagationSolution instance. We will look at the various properties and methods this class provides below, in the context of our example.

Plotting properties of the beam propagation

One of the most common use-cases is to plot a beam trace, in order to see how, in particular, the beam size and accumulated Gouy phase evolve over a path. This is as simple as calling PropagationSolution.plot() on the resulting solution. An example is given below for the call in the above section.

ps.plot(
    name_xoffsets={"ITM_lens": 10, "SRM": 10},
    name_yoffsets={"ITM_lens": 10},
);

The name_xoffsets and name_yoffsets arguments are optional and provide a way of shifting the positions of the component names on the figure (in terms of data coordinates) to avoid clashes.

Just the beam-sizes, for example, can also be plotted.

ps.plot_beamsizes(
    name_xoffsets={"ITM_lens": 10, "SRM": 10},
    name_yoffsets={"ITM_lens": 10},
);

The wavefront curvature is not plotted by default with a call to PropagationSolution.plot() but this can be plotted separately too with PropagationSolution.plot_curvatures() or along with the others via giving "all" as the first argument to this method.

Printing propagation data

As well as plotting, one can also display useful data associated with a beam propagation by simply printing the solution object. For our example here, this results in:

print(ps)

┌───────────────╥──────────┬───────────┬───────────┬───────────┬────────────┬────────────┬───────────┬─────────────────────┐
│               ║    z     │    w0     │    zr     │     w     │    RoC     │     S      │ Acc. Gouy │          q          │
╞═══════════════╬══════════╪═══════════╪═══════════╪═══════════╪════════════╪════════════╪═══════════╪═════════════════════╡
│   ITM.p1.o    ║      0 m │ 8.9093 mm │  550.22 m │ 89.907 mm │    5.58 km │  179.21 uD │        0° │ 5525.206 + 550.225j │
├───────────────╫──────────┼───────────┼───────────┼───────────┼────────────┼────────────┼───────────┼─────────────────────┤
│  ITMAR.p2.i   ║   200 mm │ 8.9093 mm │  550.22 m │  89.91 mm │  5.5802 km │  179.21 uD │   204.5u° │ 5525.406 + 550.225j │
├───────────────╫──────────┼───────────┼───────────┼───────────┼────────────┼────────────┼───────────┼─────────────────────┤
│  ITMAR.p1.o   ║   200 mm │ 8.9093 mm │  160.88 m │  89.91 mm │  1.6316 km │   612.9 uD │   204.5u° │ 1615.573 + 160.880j │
├───────────────╫──────────┼───────────┼───────────┼───────────┼────────────┼────────────┼───────────┼─────────────────────┤
│ ITM_lens.p2.i ║   200 mm │ 8.9093 mm │  160.88 m │  89.91 mm │  1.6316 km │   612.9 uD │   204.5u° │ 1615.573 + 160.880j │
├───────────────╫──────────┼───────────┼───────────┼───────────┼────────────┼────────────┼───────────┼─────────────────────┤
│ ITM_lens.p1.o ║   200 mm │ 431.39 um │ 377.18 mm │  89.91 mm │  -78.614 m │  -12.72 mD │   204.5u° │    -78.612 + 0.377j │
├───────────────╫──────────┼───────────┼───────────┼───────────┼────────────┼────────────┼───────────┼─────────────────────┤
│   ZM2.p2.i    ║   52.7 m │ 431.39 um │ 377.18 mm │ 29.868 mm │  -26.117 m │ -38.289 mD │  552.87m° │    -26.112 + 0.377j │
├───────────────╫──────────┼───────────┼───────────┼───────────┼────────────┼────────────┼───────────┼─────────────────────┤
│   ZM2.p1.o    ║   52.7 m │ 1.1751 mm │  2.7989 m │ 29.868 mm │  -71.193 m │ -14.046 mD │  552.87m° │    -71.083 + 2.799j │
├───────────────╫──────────┼───────────┼───────────┼───────────┼────────────┼────────────┼───────────┼─────────────────────┤
│   ZM1.p2.i    ║  102.7 m │ 1.1751 mm │  2.7989 m │ 8.9293 mm │  -21.454 m │ -46.611 mD │   5.8603° │    -21.083 + 2.799j │
├───────────────╫──────────┼───────────┼───────────┼───────────┼────────────┼────────────┼───────────┼─────────────────────┤
│   ZM1.p1.o    ║  102.7 m │ 6.1019 mm │  75.465 m │ 8.9293 mm │  -151.26 m │ -6.6111 mD │   5.8603° │   -80.625 + 75.465j │
├───────────────╫──────────┼───────────┼───────────┼───────────┼────────────┼────────────┼───────────┼─────────────────────┤
│    BS.p2.i    ║  172.7 m │ 6.1019 mm │  75.465 m │ 6.1621 mm │  -546.61 m │ -1.8295 mD │   44.739° │   -10.625 + 75.465j │
├───────────────╫──────────┼───────────┼───────────┼───────────┼────────────┼────────────┼───────────┼─────────────────────┤
│    BS.p4.o    ║  172.7 m │ 9.7866 mm │  281.41 m │ 9.8832 mm │ -2.0383 km │  -490.6 uD │   44.739° │  -39.622 + 281.411j │
├───────────────╫──────────┼───────────┼───────────┼───────────┼────────────┼────────────┼───────────┼─────────────────────┤
│  BSAR2.p1.i   ║ 172.77 m │ 9.7866 mm │  281.41 m │ 9.8828 mm │  -2.042 km │ -489.71 uD │   44.754° │  -39.547 + 281.411j │
├───────────────╫──────────┼───────────┼───────────┼───────────┼────────────┼────────────┼───────────┼─────────────────────┤
│  BSAR2.p3.o   ║ 172.77 m │ 6.1019 mm │  75.465 m │ 6.1618 mm │   -547.6 m │ -1.8262 mD │   44.754° │   -10.605 + 75.465j │
├───────────────╫──────────┼───────────┼───────────┼───────────┼────────────┼────────────┼───────────┼─────────────────────┤
│   SRM.p1.i    ║ 182.77 m │ 6.1019 mm │  75.465 m │ 6.1021 mm │   -9.41 km │ -106.27 uD │   52.294° │    -0.605 + 75.465j │
└───────────────╨──────────┴───────────┴───────────┴───────────┴────────────┴────────────┴───────────┴─────────────────────┘

A matrix of distances between all the optics in the propagation can also be printed using PropagationSolution.distances_matrix_table():

print(ps.distances_matrix_table())

┌──────────╥───────────┬───────────┬───────────┬───────────┬───────────┬───────────┬──────────┬──────────┐
│          ║    ITM    │   ITMAR   │ ITM_lens  │    ZM2    │    ZM1    │    BS     │  BSAR2   │   SRM    │
╞══════════╬═══════════╪═══════════╪═══════════╪═══════════╪═══════════╪═══════════╪══════════╪══════════╡
│   ITM    ║       0 m │    200 mm │    200 mm │    52.7 m │   102.7 m │   172.7 m │ 172.77 m │ 182.77 m │
├──────────╫───────────┼───────────┼───────────┼───────────┼───────────┼───────────┼──────────┼──────────┤
│  ITMAR   ║   -200 mm │       0 m │       0 m │    52.5 m │   102.5 m │   172.5 m │ 172.57 m │ 182.57 m │
├──────────╫───────────┼───────────┼───────────┼───────────┼───────────┼───────────┼──────────┼──────────┤
│ ITM_lens ║   -200 mm │       0 m │       0 m │    52.5 m │   102.5 m │   172.5 m │ 172.57 m │ 182.57 m │
├──────────╫───────────┼───────────┼───────────┼───────────┼───────────┼───────────┼──────────┼──────────┤
│   ZM2    ║   -52.7 m │   -52.5 m │   -52.5 m │       0 m │      50 m │     120 m │ 120.07 m │ 130.07 m │
├──────────╫───────────┼───────────┼───────────┼───────────┼───────────┼───────────┼──────────┼──────────┤
│   ZM1    ║  -102.7 m │  -102.5 m │  -102.5 m │     -50 m │       0 m │      70 m │ 70.075 m │ 80.075 m │
├──────────╫───────────┼───────────┼───────────┼───────────┼───────────┼───────────┼──────────┼──────────┤
│    BS    ║  -172.7 m │  -172.5 m │  -172.5 m │    -120 m │     -70 m │       0 m │ 74.78 mm │ 10.075 m │
├──────────╫───────────┼───────────┼───────────┼───────────┼───────────┼───────────┼──────────┼──────────┤
│  BSAR2   ║ -172.77 m │ -172.57 m │ -172.57 m │ -120.07 m │ -70.075 m │ -74.78 mm │      0 m │     10 m │
├──────────╫───────────┼───────────┼───────────┼───────────┼───────────┼───────────┼──────────┼──────────┤
│   SRM    ║ -182.77 m │ -182.57 m │ -182.57 m │ -130.07 m │ -80.075 m │ -10.075 m │    -10 m │      0 m │
└──────────╨───────────┴───────────┴───────────┴───────────┴───────────┴───────────┴──────────┴──────────┘

Accessing beam properties

PropagationSolution provides various methods for accessing physical properties of the beam at any point along the computed path. All of these relevant methods derive from PropagationSolution.q() (i.e. the beam parameter), so here we will explore the options available for accessing q at different locations.

The simplest case is to provide an OpticalNode instance to get the beam parameter at this node:

print(ps.q(model.BS.p2.i))

BeamParam(w0=6.1019 mm, z=-10.625 m, w=6.1621 mm, Rc=-546.61 m)

in the example above we grab the beam parameter at the input node of the second port of the beam splitter (coming from the ZM1 optic). Whilst model.BS.p2.o (i.e. the output node of this port) is technically not traced in the given path, we can still access the beam parameter at this node as the beam tracing assumes that opposite node beam parameters are the reverse of the forward node (i.e. \(-q^*\)):

print(ps.q(model.BS.p2.o))

BeamParam(w0=6.1019 mm, z=10.625 m, w=6.1621 mm, Rc=546.61 m)

One can also use the string representation of an optical node:

print(ps.q("BS.p2.i"))

BeamParam(w0=6.1019 mm, z=-10.625 m, w=6.1621 mm, Rc=-546.61 m)

Note

Any beam parameter value returned by PropagationSolution.q() will be a BeamParam instance, meaning all the various properties of this can be accessed as usual. PropagationSolution does provide some shortcuts for the key parameters (such as beam-size with PropagationSolution.w()), however, so the choice is up to the user as to whether to access like this:

ps.q(model.BS.p2.i).w

0.0061620647495966045

or like this:

ps.w(model.BS.p2.i)

0.0061620647495966045

Obtaining cumulative Gouy phases

It can be useful to see accumulated Gouy phases over segments of the traced path too, and PropagationSolution provides methods to obtain these. These phases are always given in degrees.

To obtain the total accumulated Gouy phase over the full traced path one can simply do:

ps.total_acc_gouy

52.29424393023219

Accumulated Gouy phase over a specific sequence of spaces can be retrieved with PropagationSolution.acc_gouy(), e.g:

ps.acc_gouy("lZM2_ITMlens", "lZM1_ZM2")

5.860050120196138

gives the Gouy phase accumulated from the ITM lens to ZM1 (or vice-versa). And finally, the Gouy phase accumulated up to a specific point can be obtained with PropagationSolution.acc_gouy_up_to(). This gets the cumulative Gouy phase from the starting node of the propagation up to the specified point. Similarly to accessing beam properties, this point can be an OpticalNode instance or a component or name of a component. For example, the below retrieves the accumulated Gouy phase from the ITM to BS:

ps.acc_gouy_up_to("BS")

44.739318811593286

Symbolic beam propagation

The default behaviour of propagate_beam() is to use the current values of each parameter, returning non-symbolic data in the PropagationSolution instance. Another way to use this function is to switch on the symbolic flag - resulting in symbolic expressions for all beam parameters, Gouy phases and ABCD matrices at each point in the solution.

Using the example file from above again, we can simply call:

ps_sym = model.propagate_beam(path=ITM_TO_SRM, symbolic=True)

to compute a symbolic propagated beam solution. Accessing beam properties (such as the beam parameter via PropagationSolution.q()) will now return a symbolic expression rather than just a numeric value. Note that in the case of PropagationSolution.q(), a BeamParam instance is still returned but this will be a symbolic beam parameter (as indicated by the BeamParam.symbolic flag).

Symbolic beam propagation can result in very long symbolic expressions. Rudimentary simplification routines are provided in Finesse but for large models keeping every parameter symbolic is unnecessary. If you only need specific parameters to be kept as symbols then you can specify which to keep, for example, symbolic=(‘ITMX.Rcx’,) or symbolic=(‘ITMX.Rcx’, ‘ITMYlens.f’). If simplification of symbolic beam propagation is required you must specify simplify=True, it does not happen by default. This is because simplifying when using every single parameter takes too long and is to complicated to work with.

Plotting symbolic beam propagations

Plotting of the solution is supported, with the added option of substituting parameters in the PropagationSolution.plot() call. An example is given below, where the RoC of ZM1 is changed from \(R_c = -50\) m to \(R_c = -70\) m, and the distance between the two telescope mirrors (ZM1 , ZM2) is reduced to 30 metres:

ps_sym.plot(subs={model.ZM1.Rcx: -70, model.lZM1_ZM2.L: 30});

If subs is not given, then the current value of each parameter will be used.

Animation of a symbolic PropagationSolution is also supported via the method PropagationSolution.animate(). This method expects at least one array-like parameter substitution in the subs dict in order to perform the animation over this parameter scan. For example, we can see how the beam sizes and accumulated Gouy phases change as we scan over the telescope length via:

import numpy as np

# Vary telescope length from 20 m to 70 m
tel_zs = np.linspace(20, 70, 60)
fig, axes, anim = ps_sym.animate(
    {model.lZM1_ZM2.L: tel_zs}, interval=50,
)

Once Loop Reflect

Evaluating symbolic beam properties

The true power of symbolic PropagationSolution objects comes from the ability to obtain symbolic expressions for any geometric property of the beam at any point in the traced path. Accessing these properties is performed in exactly the same way as detailed in the “Accessing beam properties” section above.

As an example, we can obtain a symbolic expression for the beam size on ZM1 with:

# Symbolic expr of w at ZM1
w_zm1_sym = ps_sym.w(model.ZM1.p2.i)

# Evaluate this with current parameter values
print(f"w_ZM1 = {w_zm1_sym.eval() / 1e-3} mm")

w_ZM1 = 8.929276136051481 mm

Similarly to the plotting routines, we can pass a subs dict to the eval method of the symbolic expression to find the beam size at ZM1 using different dependent optic parameters:

# Evaluate w_ZM1 with ZM2 RoC changed to -70 m
print(f"w_ZM1 = {w_zm1_sym.eval(subs={model.ZM2.Rcx: -70}) / 1e-3} mm")

# Evaluate w_ZM1 with ITM lens focal length changed to 100 m
# and telescope length changed to 60 m
print(f"w_ZM1 = {w_zm1_sym.eval(subs={model.ITM_lens.f: 100, model.lZM1_ZM2.L: 60}) / 1e-3} mm")

w_ZM1 = 15.378018812810746 mm
w_ZM1 = 61.29723873257198 mm

These symbolic expressions also support numpy.ndarray type arguments, allowing for evaluation over N-dimensional arrays. For example, we can evaluate the beam size at ZM1 over a range of ITM lens focal lengths and plot the results:

import matplotlib.pyplot as plt

fs = np.linspace(50, 90, 100)
w_zm1s = w_zm1_sym.eval(subs={model.ITM_lens.f: fs}) / 1e-3 # scale to mm

plt.plot(fs, w_zm1s)
plt.xlabel("ITM lens focal length [m]")
plt.ylabel("Beam size at ZM1 [mm]");

Or we could, for example, compute the Gouy phase accumulated up to ZM1 for the same parameter space:

acc_gouy_zm1_sym = ps_sym.acc_gouy_up_to("ZM1")

acc_gouys_zm1 = acc_gouy_zm1_sym.eval(subs={model.ITM_lens.f: fs})

plt.plot(fs, acc_gouys_zm1)
plt.xlabel("ITM lens focal length [m]")
plt.ylabel("Accumulated Gouy phase from\nITM to ZM1 [deg]");

This is a very powerful feature, as it leverages the speed of NumPy array calculations allowing for fast, high-dimensional grid calculations of arbitrary beam properties over any dependent model parameter.

Reverse beam propagation

Some components in a model may not have forwards and backwards couplings between ports, such as a directional beamsplitter. These elements act as a form of optical isolation, such as a faraday isolator would. This means you can only trace beams through a physically viable path, as the beam tracer only follows the connections at each component.

To get around this, it is possible to “reverse propagate” and let FINESSE know you want to get around this limitation and follow unphysical paths. Here is a simple example:

model = finesse.script.parse(
    """
l l1
dbs isolator
lens L1 f=1
m m1
link(l1.p1, 1, isolator.p1, isolator.p3, 2, L1, 3, m1)
gauss g1 l1.p1.o w0=1e-3 z=1e-2
"""
)

fwd = model.propagate_beam("l1.p1.o", "m1.p1.i")
rev = model.propagate_beam("m1.p1.o", "l1.p1.i", reverse_propagate=True)

fwd.plot()
rev.plot()

(<Figure size 576x355.968 with 2 Axes>,
 array([<Axes: ylabel='Beam size [mm]'>,
        <Axes: xlabel='Distance [m]', ylabel='Gouy phase\naccumulation [deg]'>],
       dtype=object))