Thermal lensing and deformations using Hello-Vinet

A small fraction of power is absorbed from an incident beam whenever it transmits through or reflects from a physical optical component. This absorbed power is dumped into the optic in the form of heat. Two main effects are apparent when this occurs; a thermal gradient develops within the substrate of the optic and the refractive index changes, the thermo-refractive effect; and the optic expands and deforms itself and any reflective surface, the thermo-elastic effect.

Determining how these processes affect the optical field usually requires finite element analysis tools solving the thermal diffussion and linear elastic problems for a 3D geometry. Alternatively there are also a collection of analytic solutions for on-axis cylindrically symmetric heating beams interacting with a cylindrical mirror. These formulae are known as the Hello-Vinet equations within the gravitational wave community, after the two authors who derived them. A farily comprehensive review article [32] exists which summarises the work by Vinet.

Finesse provides some of the Hello-Vinet calculations which can all be found in finesse.thermal.hello_vinet. In this module are functions for computing heating from fundamental HG00 Gaussian beams as well as more generic axisymmetric heating irradiance patterns. The former functions have a prefix _HG00.

Steady state substrate temperature

The thermo-refractive effect generates a thermal lens within the substrate of an optic. The temperature change throughout the substrate can be found using the following function:

import finesse
import numpy as np
import matplotlib.pyplot as plt
import finesse.thermal.hello_vinet as hv
from finesse.materials import FusedSilica
finesse.init_plotting()

a = 0.17 # mirror radius
h = 0.2  # mirror thickness
w = 53e-3 # spot size radius
r = np.linspace(-a, a, 50) # radial points
z = np.linspace(-h/2, h/2, 100) # longitudinal points
material = FusedSilica

T_coat_per_W, T_bulk_per_W = hv.substrate_temperatures_HG00(
    r, z, a, h, w, material
)

The returned arrays are temperature change per Watt of optical power. Therefore to get the actual temperature changes you must multiply the results by the incident power on the coating or the power going through the substrate.

plt.title("Coating heating")
plt.contourf(r, z, T_coat_per_W)
plt.xlabel("r [m]")
plt.ylabel("z [m]")
plt.colorbar(label='dT/W [K]')

<matplotlib.colorbar.Colorbar at 0x781ad608e870>

plt.title("Bulk heating")
plt.contourf(r, z, T_bulk_per_W)
plt.xlabel("r [m]")
plt.ylabel("z [m]")
plt.colorbar(label='dT/W [K]')

<matplotlib.colorbar.Colorbar at 0x781ad6102180>

Currently the code only supports HG00 heating beams.

Steady state thermal lensing

The effective thermal lens due to this temperature change throughout the substrate is computed with the integral

\[Z(r) = \frac{dn}{dT} \int T(r,z) dz\]

The above finesse.thermal.hello_vinet.substrate_temperatures_HG00() outputs can be numerically integrated, but we can also use the finesse.thermal.hello_vinet.thermal_lenses_HG00() method for a direct calculation. Using the same parameters from before we find the thermal lens with

Z_coat_per_W, Z_bulk_per_W = hv.thermal_lenses_HG00(
    r, a, h, w, material
)

plt.plot(r, Z_coat_per_W)
plt.xlabel("r [m]")
plt.ylabel("OPD [m]")

Text(0, 0.5, 'OPD [m]')

The returned 1D arrays are the optical path difference (OPD) in meters. This is the extra length difference the beam will experience on a single pass through the mirror. Again this has not been scaled for any incident power yet.

Steady state thermal displacements

Displacements within an optic are also caused by the thermal expansion of a material as it heats up. An Incident HG00 beam deform the surfaces of the mirrors as well as slightly expand the substrate.

Substrate displacements from coating absorption can be computed with:

w = 53e-3
a = 170e-3 # radius
h = 0.2 # thickness
z = np.linspace(-h/2, h/2, 1000)
r = np.linspace(-a, a, 100)
# z displacement throughout the substrate
U_z_coat = hv.substrate_thermal_expansion_depth_HG00(r, z, a, h, w, FusedSilica)

plt.contourf(r, z, U_z_coat)
plt.colorbar(label='z displacement [m]')
plt.xlabel("r [m]")
plt.ylabel("z [m]")
plt.title("Coating thermal deformtion")

Text(0.5, 1.0, 'Coating thermal deformtion')

We can also just plot the surface deformation on the absorbing side.

w = 53e-3
a = 170e-3 # radius
h = 0.2 # thickness
r = np.linspace(-a, a, 100)

U_s_coat = hv.surface_deformation_coating_heating_HG00(r, a, h, w, FusedSilica)

plt.plot(r, U_s_coat)
plt.ylabel('Surface displacement [m/W]')
plt.xlabel("r [m]")
plt.title("Coating absorption surface deformation")

Text(0.5, 1.0, 'Coating absorption surface deformation')

Care should be taken when using these deformation calculations in simulations to ensure that the correct coordinate system is being used. For example, ITM and ETM when modelled with mirrors in Finesse, typically the ITM HR surface is on the port 2 side, whereas for the ETM it is on the port 1 side. Positive z direction for mirrors are for the surface normal on the port 1 side.

In such cases the above hv.surface_deformation_coating_heating_HG00() should give the correct displacement for a beam incident on the port 2 side of a mirror, making it the correct displacement for the ITM mirror.

Axisymmetric heating profiles

The functions discussed above also have generic versions which compute the thermal lensing and deformation for generic axisymmetric heating profiles. The first step is to compute the Fourier-Bessel expansion of the radial heating profile. With this data you can then pass it directly to the Hello-Vinet functions to compute the resulting effects.

Here is an example computing how a more complex heating pattern is decomposed:

import finesse.materials
import finesse.thermal.hello_vinet as hv
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import eval_hermite

finesse.init_plotting()
material = finesse.materials.FusedSilica
a = 0.17
h = 0.2
w = 53e-3
r = np.linspace(0, a, 101)

# A non-normalised 5th order hermite radial distribution
E = eval_hermite(5, np.sqrt(2)*r/w) * np.exp(-(r/w)**2)
I = E*E
plt.plot(r, I, label='I(r)')

# perform Fourier-Bessel decomposition
data = hv.get_p_n_s_numerical(I, a, 10, material)
plt.plot(r, hv.eval_p_n_s_numerical(data), ls='--', lw=2, label='s_max=10 fit')

# perform Fourier-Bessel decomposition with
# a higher order decomposition
data = hv.get_p_n_s_numerical(I, a, 20, material)
plt.plot(r, hv.eval_p_n_s_numerical(data), ls='--', lw=2, label='s_max=20 fit')

plt.legend()
plt.xlabel("r [m]")
plt.ylabel("Intensity [Wm^-2]")
plt.title("Fourier-Bessel decomposition of irradiance")

Text(0.5, 1.0, 'Fourier-Bessel decomposition of irradiance')

It is important to check that your heating profile is adequately described by the Fourier-Bessel decomposition. Higher spatial variance requires a larger number of roots to be found (the s_max parameter). When using low spatial sampling it may also be necessary to experiment with the Newton-Cotes weighting order for more accurate numerical integration. Higher spatial frequency changes in the heating profile will not decompose into bessel functions particularly well. This is not so much of an issue as naturally heating profiles tend to be smooth and slowly varying.

Next we can use this decomposition data in the generic forms of the thermal equations. For example, computing the thermal lensing in a substrate:

W_coat, W_bulk = hv.thermal_lenses(
    data, h
)
plt.plot(r, W_coat/1e-6, label='Coating heating')
plt.plot(r, W_bulk/1e-6, label='Substrate heating')
plt.xlabel('r [m]')
plt.ylabel('OPD [um/W]')
plt.legend()

<matplotlib.legend.Legend at 0x781ad242f680>

or computing the temperature distribution through the substrate:

z = np.linspace(-h/2, h/2, 300)
T_coat, T_bulk = hv.substrate_temperatures(
    data, z, h
)
plt.figure()
plt.contourf(r, z, T_coat, levels=20)
plt.colorbar(label='Substrate temperature [K/W]')
plt.xlabel("r [m]")
plt.ylabel("z [m]")
plt.title("Coating heating")

plt.figure()
plt.contourf(r, z, T_bulk, levels=20)
plt.colorbar(label='Substrate temperature [K/W]')
plt.xlabel("r [m]")
plt.ylabel("z [m]")
plt.title("Substrate heating")

Text(0.5, 1.0, 'Substrate heating')