Global optimization of multilayer dielectric coatings for precision measurements

We describe the design of optimized multilayer dielectric coatings for precision laser interferometry. By setting up an appropriate cost function and then using a global optimizer to find a minimum in the parameter space, we were able to realize coating designs that meet the design requirements for spectral reflectivity, thermal noise, absorption, and tolerances to coating fabrication errors.


Introduction
Interferometric gravitational wave detectors such as Advanced LIGO require mirrors with dielectric coatings that satisfy multiple requirements on reflectivity at one or more wavelengths, surface electric field, absorption within the coating, and thermal noise [1].For the next generation of detectors, it is anticipated that there will be tighter requirements on these specifications.Furthermore, it is desirable that the coating design chosen will have minimal sensitivity to manufacturing tolerances.
In this paper, we describe the construction of a cost function that quantifies how closely a given coating design satisfies the multiple requirements on it.Once this cost function has been constructed, it can be minimized with a numerical optimization algorithm, such as Particle Swarm Optimization [2] or Differential Evolution [3,4], to compute the coating design (i.e. the set of thicknesses of various layers in a dielectric stack of alternating high-and low-index materials) that gives the (global) minimum value in the allowed parameter space.Since most numerical optimization algorithms are designed to handle a scalar cost function, we convert our vector (multi-objective) cost function to a scalar number by taking the scalar product of it with a weight vector.The latter allows us to quantify which design objectives are more important than others.
Once an optimal solution has been arrived at, we verify its sensitivity to small perturbations in layer thicknesses as well as various assumed model parameters using Monte-Carlo (MC) analysis.This approach allows us to evaluate the relative performance of coating designs that vary along 'hyper-parameter' axes, such as the choice of number of layer pairs that make up the coating.This approach of mapping a complex design problem into a scalar objective function minimization problem can be readily generalized to several other problems in the field of Gravitational Wave Astronomy.It is also of interest to the broader category of experiments that use optical cavities for precision measurements in the atomic and molecular optics community.

Requirements on multilayer dielectric coatings
In this section, the requirements on dielectric coatings of the type considered in this paper are described.Although the method is easily generalized to a multilayer stack with arbitrary arXiv:2110.13437v2[physics.optics]12 Dec 2023 refractive index profile, we only consider cases in which the stack is comprised of alternating layers of two dielectrics, with refractive indices  1 (low-index) and  2 (high-index).The coating is deposited on a substrate with refractive index  sub .These considerations are typical in the field of laser interferometry, but are also of relevance in various other fields such as nanoscale optomechanics.Strategies to increase the numerical efficiency of computing these properties are also briefly discussed.

Spectral reflectivity
The primary requirement on a dielectric coating design is the power reflectivity at wavelengths of interest, ().For a given dielectric coating stack with  interfaces, this may be calculated in several ways.We opt to do the computation of the amplitude reflectivity, Γ recursively, using the relation [5] where the index  = ,  − 1, ..., 1, with  = 1 corresponding to the interface of the coating with the incident electromagnetic field.In Equation (1), where   is the physical thickness of,   is the (wavelength dependent) refractive index of, and   is the angle of incidence into the  − th layer.The recursion relation Equation ( 1) is initialized with , where   sub is defined by Equation (2b) for the substrate onto which the dielectric stack is deposited.The power reflectivity of the stack may then be computed as  = |Γ 1 | 2 .A typical coating design requirement will specify (, , polarization).

Surface electric field and absorption
In optical cavities in which the circulating power is high, absorption in the dielectric layers becomes important for a number of reasons.The coatings have to be able to withstand the absorptioninduced thermal heating for the highest expected incident electromagnetic field intensity.In applications, where the cavity mirrors have to be maintained at cryogenic temperatures, the requirement becomes even more stringent as the rate at which heat can be extracted from the optic will set the maximum permissible absorption in the coating [6].
Absorption is quoted as a dimensionless fraction of the incident power which is converted to thermal energy in the coating.We evaluate absorption by first determining the square of the electric field as a function of penetration depth, z, normalized by the incident electric field, [7].The total absorption,  T in a coating of thickness  is with (z) describing the bulk absorption of the materials used in the coating, for which measured values are available.In practise, the integral in Equation (3) can be evaluated numerically.However, this is a computationally expensive operation.A good proxy for use in an optimizer is the value of the electric field transmitted into the first layer of coating, given by [8] | ì ì  +  is set by the incident power density, which varies in different applications.Since Γ 1 has to be computed for evaluating the power reflectivity, the additional computational overhead is minimized.In order to minimize ì  surface , and hence the absorption, it is conventional to add a nearly half-wavelength optical thickness 'cap' of dielectric material to a high-reflectivity (HR) coating.

Thermo-optic and Brownian noise
Sensitivity of the current generation of laser interferometric gravitaitonal wave detectors is expected to be limited by the Brownian thermal noise of the coatings used.These consist of up to 20 pairs of alternating layers of SiO 2 (low-index material) and Ta 2 O 5 (high-index material), or 40 -50 pairs in materials with low index contrast (e.g. the pair of crystalline materials AlGaAs/GaAs [9]).In order to improve the sensitivity of future generation of detectors, an active area of research pursued in the last decade is the development of alternative dielectrics with which coatings that can meet the power reflectivity requirements can be developed.
Another application in which the Brownian noise of dielectric coatings can be a limiting noise source is the development of ultrastable frequency reference cavities [10].The frequency stability of lasers used in precision metrology experiments are often stabilized to such reference cavities, and efforts are underway to identify alternative dielectrics so that even more stable reference cavities can be constructed.
There is also a second effect which has to be simultaneously considered, the 'thermo-optic' noise, which arises as a result of the temperature dependence of the refractive index of the dielectrics, and their thermally driven length fluctuations.With clever design, it is possible to suppress this noise contribution by coherent cancellation of the two effects [11].

Immunity to small perturbations in assumed model parameters
In evaluating the coating properties, assumptions are made about the layer thicknesses, refractive indices, dispersion, bulk absorption, and mechanical loss angle of the dielectrics.Additionally, uncertainties may exist, for example, in the assumed value of the angle of incidence.Since parameters of interest such as the power reflectivity of the coating are functions of these parameters, any uncertainty in them (due to manufacturing process limitations, measurement errors etc.) propagate through to the performance of the manufactured coating.In order to meet the tight tolerances on these parameters, it is desirable to choose (for fabrication) the coating design whose sensitivity to small errors in these model parameters is low.
We address this using a two step approach.First, during the optimization stage, the (numerical) derivative of coating properties with respect to model parameters are used in constructing the cost-function to be minimized.Secondly, the sensitivity of a given design to small (≈ 1%) perturbations to the assumed model parameters is evaluated explicitly.The errors themselves are assumed to be i.i.d., and are drawn from an uncorrelated multivariate zero-mean Gaussian distribution whose standard deviation is chosen to be 0.5% of the optimized value of the parameter.For example, errors in the thickness of the -th layer,   , is sampled from the distribution ( opt  and  opt  being the value of the thickness of the -th layer that best achieves the design objectives.The emcee package [12] is used to generate ≈ 1 × 10 5 samples, and the corner package [13] provides a convenient way to visualize the results.A confidence interval on meeting specifications within tolerance can thus be stated.

Adjust coating layer thicknesses
No Global minimum or termination criteria reached?Fig. 1: Overview of the workflow adopted.As shown in Figure 1a, the process starts by identifying the requirements, and encoding these in a cost function.The cost function is then minimized using a global optimizer, such as MATLAB's particle swarm optimization toolbox [2] or SciPy's Differential Evolution [3,4].In order to choose between multiple candidate solutions, Monte-Carlo sensitivity analysis is used to choose the solution with the least sensitivity to model parameters.Figure 1b shows the advantage of defining the cost function in a piecewise manner.

Numerical optimization techniques for coating design
The methodology adopted for optimizing the dielectric coating design is schematically illustrated in Figure 1a.

Cost function construction and minimization
The parameter space over which the optimization is done is the set of coating layer thicknesses, {  }, and in the problems we consider, can be O (100)-dimensional.Furthermore, the values   can take are constrained to not be too small (for practical manufacturing reasons) and not larger than the longest optical wavelength of interest, .Once a set of requirements,   , and weights that reflect their relative importance,   have been arrived at, these are encoded into a cost function, C(  ,   ,   ).The weights are necessary to convert the multiple objectives of the optimization problem into a scalar which can be minimized using a numerical optimzation algorithm.The mathematical statement of the optimization problem amounts to min We chose MATLAB's Particle Swarm Optimization (PSO) toolbox for implementing Equation (5).Typical runtimes on a machine that can perform 45 GFlops is O (10) minutes.Subsequently, we also found similar performance could be realized using SciPy's Differential Evolution optimizer [3], which has the added advantage of being an open-source utility.The functional form of the individual terms contributing to the cost function,   , was deliberately constructed in a normalized, piecewise manner.Normalization was necessary in order to compare costs from different requirements.Furthermore, since PSO looks for a globally optimal solution, it is likely that individual particles will traverse regions of high cost in the multi-dimensional parameter space.Defining the cost in a piecewise manner preserves a derivative term that allows the PSO functions to converge to a global minimum, but varied quadratically near the desired values and only logarithmically above some threshold value.As illustrated in Figure 1b, this prevents the cost function from blowing up to large values for poor candidate solutions, as compared to other ways of specifying the error such as the L2− and L1− norms.Other functional forms, such as sinh −1 |  −     | can also be used to achieve the same goal.

Monte-Carlo sensitivity analysis
As there are "hyper-parameters" that do not directly enter the cost function, such as the number of layers in a coating, it is possible to have the optimization algorithm yield multiple candidate solutions.We wish to choose the simplest solution (from a fabrication standpoint) that is least sensitive to small perturbations in model parameters, and meets the coating requirements within specified tolerances.For the cases considered in this paper, the perturbations used are summarized in Table 1.

Application to the inverse problem
An interesting application of the Monte-Carlo approach is to apply it to the inverse problem of inferring the optical thicknesses of a dielectric coating, given a noisy measurement of its spectral reflectivity as a function of the wavelength, .The problem becomes computationally expensive to evaluate if the dimensionality is unrestricted -that is, if we allow the physical thicknesses and refractive indices of all layers to be arbitrary.However, in practise, the dimensionality of the problem can be reduced.For concreteness, consider an optical coating composed of two dielectrics, SiO 2 and Ta 2 O 5 .Assume the coating is built up with 19 repeated identical bilayer pairs, with the top bilayer pair having a different thickness for reducing the surface electric field amplitude.Furthermore, the dispersion of the dielectrics composing the coating are well characterized, and so may be taken as fixed.In this example, the task then amounts to the following -given the power transmissivity  (), can we infer 4 numbers:  1 , the thickness of the top layer of SiO 2 ,  2 , the thickness of the next layer of Ta 2 O 5 , and [ 3 ,  4 ], the thicknesses of the repeated bilayer pair?Applying this approach to the Harmonic Separator described in Section 4.1, we infer thicknesses for the constructed coating that are within manufacturing tolerances.The modelled spectral reflectivity curve for the inferred coating is in good agreement with the measurement, as shown in Figure 2a.
Having validated the technique, we applied it to a more complex problem, namely the Advanced LIGO ETM, results for which are shown in Figure 2b.At shorter wavelengths, we found some deviation between the model constructed using inferred layer thicknesses, and the measured datathis was due to the unavailability of an accurate dispersion model for the dielectrics.Nevertheless, the technique yielded excellent results at the primary wavelength of interest, 1064 nm.Table 1: Summary of perturbations (assumed to be uncorrelated, and hence, sampled from a multi-dimensional Gaussian probability distribution with diagonal covariance matrix).These were motivated by considering realistic uncertainties on these parameters.

Model parameter Uncertainty [%]
Physical thickness of layers, Δ  † 0.5 Refractive indices of layers, Δ 1 and Δ 2 0.5 Angle of incidence, Δ 1 ‡ 1 † Due to the nature of the manufacturing process, all layers are expected to have identical fractional uncertainties in physical thickness.‡ Applicable only for non-normal incidence cases.

Case studies and results
Table 2 summarizes the requirements for a few case studies to which this methodology was applied.More details about the individual design requirements and results from the optimization runs are presented in the following subsections.---* AlGaAs is used to collectively refer to alternating layers of Al 0.92 Ga 0.08 As (low-index material) and GaAs (high-index material).†  1 = 1064nm,  2 = 532nm.‡ For noise requirements, numbers quoted are amplitude spectral densities at 100 Hz. ∧ Both polarizations are degenerate for normal incidence, as is the case for the ETMs.**  1 = 2050 m, with -Si:SiN bilayers on crystalline Silicon substrate at 124 K

Harmonic separator
The aLIGO interferometers use multiple wavelengths of laser light to sense and control interferometric degrees of freedom of the suspended optical cavities [16,17].In this case, the objective was to design a harmonic separator that allowed extraction of light at the second harmonic,  2 = 532nm, from a folded optical cavity, while preserving high reflectivity for the fundamental light field at  1 = 1064nm.Furthermore, since the expected angle of incidence on this optic was ≈ 41.1  , the design had to meet the  and  specifications for both s-and p-polarizaitons at 532nm, while only p-polarization was of interest at 1064nm.
Figure 3a shows the spectral reflectivity of the optimized coating design.Figure 3b compares the measured performance of a harmonic separator fabricated with layer thicknesses generated using this optimization routine.The measured spectral properties meet the design requirements, and are consistent with the tolerance analysis presented in Figure 3c, which evaluates the sensitivity of the design to small perturbations, as described in Section 3.2.(c) Fig. 3: Performance of the optimized harmonic separator coating design.Figure 3a indicates the wavelengths of interest,  1 (orange) and  2 (green), with dashed vertical lines.Figure 3b shows the measured performance of a harmonic separator fabricated with layer thicknesses generated using the optimization routine described in the text.Figure 3c shows the robustness of the design choice of 20 layer pairs in the dielectric stack to small variations in assumed model parameters.

HR cavity mirror coating for a gravitational wave detector
The aLIGO interferometers are designed to have Fabry-Pérot arm cavities with finesse ≈ 450, for which the output coupling mirrors of the Fabry-Pérot arm cavities (referred to as the End Test Mass, ETM) is required to have  = 5 ± 1ppm at  1 = 1064nm [18].It is also required to have  = 96.8% at  2 = 532nm to facilitate sensing and control of the arm cavity length using an auxiliary laser wavelength during the lock-acquisiton process [16,17].Additionally, since thermally driven microscopic fluctuations in the coating's optical and physical thickness is expected to limit the sensitivity of the instrument in the 100Hz-1kHz frequency band, the chosen coating design's thermo-optic (TO) and Brownian noise spectral densities should not exceed specified thresholds in this frequency range.For SiO 2 /Ta 2 O 5 coatings of the type used on the aLIGO optics, the mechanical loss angle,  dielectric is ≈ 17× larger for Ta 2 O 5 than for SiO 2 [19], and so the total thickness of the latter in a given coating design dominates the thermal noise contribution.Hence, the design should be optimized for minimum thermal noise, while still meeting other requirements.
Finally, the circulating power in these Fabry-Pérot cavities is expected to be O (1MW) during high power operation.The coating should be designed with a safety factor such that it is not damaged under these conditions.One possible damage mechanism is that residual particulate matter on the optic's surface gets burnt into the coating.In order to protect against this, the coating has to be optimized to have minimum surface electric field.
With these requirements as inputs to the optimization problem, we ran the particle swarm and sensitivity analysis and obtained a set of layer thicknesses.Figure 4 shows the performance of Fig. 4: Performance of the optimized aLIGO ETM coating design.Here, we compare the robustness of the optimized design (red) and the as-built aLIGO ETM (blue), to small variations in assumed model parameters.The superior performance of our optimized coating, with respect to the design goals, is evident.
the optimized coating, and compares it to the as-built aLIGO ETM.The superior performance of our optimized design is evident -in particular, the transmissivity at 532 nm is much more robust in our design1.We assume  Ti:Ta 2 O 5 = 3.6 × 10 −4 -this value is continually being revised by better measurements [19], but since the same value is used to compare our design to the as-built aLIGO ETM, the conclusion that our design performs better remains valid.

HR coatings for next-generation detectors
As mentioned in Section 2.3, alternative dielectric materials to SiO 2 /Ta 2 O 5 are being considered in an effort to reduce the coating Brownian noise and improve the sensitivity of next-generation laser interferometric GW detectors.A promising alternative for wavelengths around 1 m are alternating layers of low-index Al 0.92 Ga 0.08 As and high-index GaAs crystalline dielectrics shown to yield up to a 4-fold detector sensitivity improvement [20].On the other hand, cryogenic Silicon detectors expected to operate at longer wavelengths [21] have recently considered low-index Si 3 N 4 (silicon nitride) and high-index -Si (amorphous Silicon) to improve the coating thermal noise limit by over an order of magnitude relative to SiO 2 /Ta 2 O 5 at 124 K [14,22,23].
While the coating Brownian noise may be significantly lower due to the low mechanical loss in these thin films, the overall thermal noise has to take into account both the Brownian 1In Advanced LIGO's first observing run, lock acquisition was made more difficult by the fact that the reflectivity of the ETM at 532 nm was ≈ 65 %, while the design goal was 97%. Figure 4 shows that  532 nm for the as-built ETM coating design is indeed very sensitive to small errors in the model parameters.Fig. 5: Figure 5a and Figure 5b show the sensitivity analysis of the optimized AlGaAs:GaAs and Si:SiN coatings performance respectively.Below, Figure 5c and Figure 5d show the normalized square magnitude of electric field |E(z)/E +  | 2 inside the AlGaAs:GaAs and Si:SiN coatings respectively, along with the physical thickness profile of the individual dielectric layers after optimization.Dashed vertical lines indicate boundaries between low-and high-index materials.and thermo-optic noise contributions.The latter can be minimized by coherent cancellation of thermorefractive and thermoelastic effects.This is done by including a penalty for the TO noise at a representative frequency (we choose 100 Hz) in the cost function.While the absorption is not explicitly included in the cost function that is minimized by PSO, we include it in the MC sensitivity analysis, and confirm that the likelihood of it lying within the acceptable range of ≤ 1ppm is high, even if there are small deviations in assumed model parameters.The overall performance of the optimized crystalline coating is shown in Figure 5a, while a similar sensitivity plot is shown in Figure 5b for the optimized cryogenic Si:SiN coating.Figure 5c and Figure 5d show the variation of the electric field inside the crystalline and cryogenic coatings respectively, along with their layer thickness profile.

Survey of past work, and future outlook
There have been many proposed optimization algorithms for realizing dielectric coating thicknesses [24][25][26][27].The general algorithm is to define a cost function based on the design objectives, and then use some algorithm to find the global minimum value of this cost function in the allowed parameter space.Popular choices for global cost function minimization are the genetic algorithm, needle optimization, and swarm optimization algorithms [28].Compared to these past works, our approach can easily handle these multiple design goals in a weighted manner.Moreover, our algorithm preferentially selects coating designs that are minimally sensitive to manufacturing tolerances, thereby decreasing the probability of fabricating a coating that does not meet the design specifications.Our optimization code has yielded comparable or superior coating designs relative to commercial coating design software, and since it is built using Python, does not incur expensive licensing fees.Development of a GUI will make the experience even smoother for a casual user, while more advanced problems of interest in the coating community can be set up easily given the modular nature of the algorithm.
During the writing of this manuscript, others in the thin-film community have developed complementary tools that can readily be integrated with the work presented here.In particular, [29] has leveraged recent developments from the machine-learning community to dramatically reduce the computational time for evaluating the recursive relations described in Section 2.1.Adapting that software suite to the inverse problem described in Section 3.2.1 could offer accurate inference of thin-film structures from spectral measurements for even more complex coatings than the Advanced LIGO ETM.Designing anti-reflective (AR) coatings for transmitting carefully prepared quantum states of light with minimal loss in optical systems is another case-study we are currently investigating.

Conclusion
We have presented a method to find a globally optimum solution to multilayer dielectric coating design problems such that multiple competing objectives are satisfied.We only use standard libraries available with MATLAB's PSO toolbox and Python for the implementation of this algorithm.PSO was chosen as the optimization algorithm given its success in similar problems.We found that using different optimization algorithms, such as Simulated Annealing, did not result in significant improvements in computational time.The performance of our optimized coatings are superior to conventional designs such as the "quarter-wave stack", and are robust to small perturbations in assumed model parameters.We have demonstrated the successful application of this technique to multiple design problems with varying requirements, including a harmonic separator, an HR coatings for aLIGO, and next-generation detectors using crystalline or cryogenic thin films.While the numerical efficiency of the algorithm may be improved, in its current iteration, our code runs in a reasonable amount of time on a modern multi-core machine capable of ≈ 50 GFlops.We are presently in the process of modifying the code to use Python for the global minimization step as well, so that the entire design process may be done using open-source software.

Fig. 2 :
Fig. 2: Inferring the coating structure from a spectral reflectivity measurement for the harmonic separator described in Section 4.1 (left), and the Advanced LIGO ETM (right).In both cases, critical wavelengths for the design are indicated by dashed vertical lines.

E
Surface [V/m]