Constraining temperature distribution inside LIGO test masses from frequencies of their vibrational modes

Thermal distortion of test masses, as well as thermal drift of their vibrational mode frequencies, present a major challenge for operation of the Advanced LIGO and Advanced VIRGO interferometers, reducing optical efficiency, which limits sensitivity and potentially causing instabilities which reduce duty-cycle. In this paper, we demonstrate that test-mass vibrational mode frequency data can be used to overcome some of these difficulties. First, we derive a general expression for the change in a mode frequency as a function of temperature distribution inside the test mass. Then we show how the mode frequency dependence on temperature distribution can be used to identify the wavefunction of observed vibrational modes. We then show how monitoring the frequencies of multiple vibrational modes allows the temperature distribution inside the test mass to be strongly constrained. Finally, we demonstrate using simulations, the potential to improve the thermal model of the test mass, providing independent and improved estimates of important parameters such as the coating absorption coefficient and the location of point absorbers.


I. INTRODUCTION
During Advanced LIGO's [1] first and second observing runs, about 100 kW of optical laser power was circulated in the Fabry Pérot arm cavities the interferometers [2]. During Observation Run 3, 200-250kW was circulating in the arm cavities [3]. It is planned that this power will increase to 0.5−1MW [4]. The heating of mirror surfaces of the test masses associated with this circulating power presents a significant technical challenge, since the thermal deformation of mirror surfaces leads to the loss of optical efficiency. Optical efficiency is reduced by increased scattered light losses from non-uniform absorption on the mirrored surfaces thermally deforming the surface and by reduced optical coupling between cavities as the beam is altered by thermal lensing. The reduced optical efficiency ultimately leads to a loss of interferometer sensitivity [5].
To address this problem, extra heating is applied to the test masses using specially positioned ring heaters and compensation plates. This is done in such a way that the thermal distortion caused by the ring heaters and compensation plates partially compensates that caused by the laser beam [6,7].
The test masses support a large and complex spectrum of vibrational modes in the frequency range 5-100kHz. The frequency of these modes depends on the temperature distribution inside the test mass. Some of these modes are the drivers of parametric instability [8,9], the control of which was limited by thermal transients [10,11]. Therefore it is useful to monitor the 3-dimensional temperature field inside each test mass for optical efficiency and parametric instability control.
In this paper we show that components of this temperature field can be measured in real time by monitoring the frequencies of multiple vibrational modes of the test masses. Some effort has already been spent designing and implementing a system that monitors small changes in mode frequencies [12,13]. It was shown that hundreds of vibrational modes are visible at the interferometer output at quiescent amplitudes.
These measurements can by extension allow estimates of the thermal distortion of the test mass mirror surfaces and distortion in thermo-optic lens in transmission of the test mass. Hartmann wavefront sensor [14] are currently used to monitor wavefront distortion in the test masses. The method proposed here compliments wavefront sensors, providing independent parameter estimates and information from the temperature field dimension along the optic axis.
The plan of the paper is as follows. In the next section we develop the mathematical formalism for computing the changes in mode frequencies. In section III the formalism is tested against a COMSOL [15] eigen-frequency analysis. In section IV we show how to use the frequency changes to make inferences about the temperature distribution inside the test masses and we discuss the limitations of these inferences due to symmetries of the test masses. In section V the estimated temperature distribution from a realistic scenario is examined and in sec-tion VI a Bayesian method for refining test mass thermal model parameters is described.

A. The preamble: linearity
The changes in the mode frequencies δω i are linear functions of the changes in the temperature inside the test mass, δT (r) (here and onwards bold-faced letters denote three-dimensional vectors). Mathematically this can be expressed as follows: where ρ(r) is the density, and functions f i (r) are form factors that will be discussed in the next subsection. It is convenient to introduce an inner product between functions, and similarly between vector fields: The factor ρ(r) ensures that the integral is restricted to the test mass volume, and as will be seen below, is useful for expressing orthogonality relations between the test mass mode displacements. Equation (2.1) can be written simply as .

B. Computation of the formfactors fi(r)
Consider a vector Langrangian displacement ξ(r, t) of the test mass from its position of rest. In the linear approximation, the elasto-dynamic equations of motion are whereL is the operator representing the elastic restoring force and given bŷ where ǫ kl = (ξ k,l + ξ l,k )/2 (2.7) is the shear tensor, c mnkl is the elasticity tensor, and is the elastic stress tensor. Here the Einstein convention of summing over the repeating tensorial indices is assumed. A normal mode with angular frequency ω i is characterized by the wavefunction ξ (i) (r) that satisfies the following eigenequation: Importantly, the normal modes satisfy orthogonality relation Consider now a perturbation δL due to the change in temperature of the test mass: Here we take into account the fact that the elasticity tensor is temperature-dependent. Strictly speaking, there is another contribution to the change inL that is due to thermal expansion of the test mass. However, the thermal expansion coefficient LIGO test mass substrate is very small. Numerically, it is about 0.003 × d log E/dT , where E is the Young modulus of the fused silica test mass substrate. Therefore we safely neglect the thermal expansion effect as subdominant.
Consider now the first order perturbation theory for Eq. (2.9), with the perturbed elasticity operatorL + δL and perturbed normal-mode wavefunctions ξ (i) + δξ (i) . This gives (2.12) By requiring that the perturbed eigenfunction has the same norm as the unperturbed one, we impose an extra constraint We now multiply Eq. (2.12) by ξ (i) , and integrate over the volume. Using the orthogonality relations Eqs (2.10) and (2.13), and the self-adjointness of [1/ρ(r)]L, we get (2.14) Therefore, the change of the mode's angular frequency is given by We are now ready to determine the formfactor f i (r). To acheive this, we write down the numerator of the above equation explicitly as an integral over volume: The last step is obtained by integrating by parts, using Gauss' theorem, recalling that σ mn = δσ mn = 0 at the surface of the test mass, and using the symmetry of the elasticity tensor with respect to the indices m and n (the latter insures that the stress tensor is symmetric). From this expression, we conclude that the formfactor is given by where the normalization factor is given by where E i is the total energy of the mode. It is worth noting that c mnkl ǫ mn ǫ kl = 2U (r) (2.18) where U (r) is the energy density of elastic deformation.
For an isotropic medium such as fused silica glass, where Y is the Young modulus, µ is the shear modulus, and ǫ s ik is the incompressible part of the shear, A simple way of rewriting the formfactor in Eq. (2.16), that may be handy in the context of using materials engineering packages like COMSOL or ANSYS, is as follows: where the notation implies that the partial derivative with respect to temperature is evaluated with the mode displacement ξ (i) (r) being held constant. This completes our computation of the formfactors.

III. NUMERICAL TEST AND A PROPOSAL FOR PRACTICAL MODE IDENTIFICATION
To validate the form factor solution of Eq. (2.21) it is compared to a finite element model eigen-frequency analysis performed with COMSOL [15]. The model used is  Table I and the model geometry is displayed in Figure 1 The form factors are calculated based on the strain distribution U (i) and total energy E i of a COMSOL eigenfrequency analysis of the test mass in the ambient (291K) temperature thermal state. Form factors and the total mode displacement are shown in Figure 15 in Appendix A for a selection of modes.
An analytically described change in temperature distribution is defined for the purpose of this test. The change in temperature distribution is defined by Zernike polynominal Z 3 1 across the circular surface of the test mass and a uniform distribution through the depth (z) of the test mass. Two eigen-frequency analyses are run in COMSOL, one at an ambient temperature of 291K and the second with the additional change in temperature distribution to produce a two sets of eigen-frequencies.
In conjunction, Eqs (2.4) and (2.21) are used to calculate the expected change in eigen-frequency for that same temperature distribution. The results in Figure 2 show very good agreement between the analytical expression and the COMSOL simulation. The identification of the mode shape of observed resonances presents a challenge. Parametric instabilities [8] of vibrational modes with frequencies as high as 47.5 kHz have been observed at LIGO [13]. At this frequency the mode density is high, resulting in several candidate modes that could potentially be causing the instability. Furthermore, theoretical calculations show that with increased circulating power there might be instabilities caused by modes with frequencies as high as 90 kHz [16]. (The recent installation of acoustic mode dampers [17] makes high frequency instabilities a lot less likely.) Knowledge of the mode shape is required to design active control schemes that apply forces to the test masses [10] or optical feedback [18] [19]. Currently modes are identified by comparing observed resonant frequencies with those computed using finite element modelling. Confident mode identification is currently limited to 17 kHz. At higher frequencies, imprecise knowledge of the elasticity parameters of fused silica produce large enough errors such that confusion between modes is a serious issue. The analytical expression for change of mode frequencies as a function of temperature presented here presents a new tool for mode identification. This could work as follows: 1. a well-controlled thermal transient perturbation is applied to the test mass, and the internal temperature distribution is computed as a function of time using finite element modelling. 2. The transient change in mode frequencies can be calculated as a function of time using the formalism presented above or finite element modelling. 3. These are compared and matched to the measured transient frequency changes in the monitored modes. By following this procedure, mode identifications can be confirmed. A simulated example of such a confirmation is shown in Figure 3 where 3 modes (colored green) have been deliberately misidentified by switching their indices. Simulated mode frequencies on the vertical axis are compared with mode frequencies calculated with the analytical expression of Eq. (2.1). In this case the temperature field is simulated in COMSOL as the steady state for 1W applied ring heater power. The COMSOL simulation includes thermal expansion and a uniform 0.1 mHz measurement noise has been added to the COMSOL simulated eigen-frequencies. The three misidentified modes can be clearly identified as outliers. The correct mode identification is critical for active control of parametric instability and is also required to make inferences about the temperature distribution from measurements of the eigen-frequencies of the test mass.

IV. CONSTRAINING THE TEMPERATURE FIELD INSIDE THE TEST MASS
One might suppose that if one is able to measure the temperature-induced frequency shifts of all of the vibrational modes to arbitrary precision, one should be able to reconstruct the 3-dimensional temperature perturbation inside the test mass. This would be an unprecedented fit for experiments with solids as far as we know. However, as we explain below, this strategy runs into problems because the form-factors f i (r) do not necessarily form a complete basis for all of the possible temperature perturbations; we show this explicitly for the case when the test mass has a reflection symmetry, as they in fact, do. We begin however in the next subsection by considering the conceptually simple case where the formfactors do form a complete basis and the temperature perturbation can, in principle, be measured.
A. Case of fi(r) forming a complete basis . Completeness allows us expand δT (r) in a series: Here we use the Einstein convention, where the summation of repeated indices is assumed. In general, one expects the functions f i (r) to be linearly independent, however in some cases where a high degree of symmetry is present, it may turn out that this is not so. In such situation, one needs to restrict the series above to a linearly independent subset of functions spanning the whole function space, so that the expansion is unique. Substituting the expansion above into Eq. (2.4) results in a matrix equation One therefore has where C −1 ij are the elements of C −1 . Since one monitors only finite amount N of the normal modes, in practice one should restrict C ij to be the N -dimensional square matrix.
B. Incompleteness of fi(r) due to symmetry of the test mass We do not in fact have a mathematical proof that the set of functions f i (r) is ever complete for a generic shape of the test mass, although intuitively it seems likely. However, a practically important counterexample is the case when the test mass has a reflection symmetry, say z → −z. In this case the vibrational modes have either odd or even parity with respect to z, but because it is the elastic energy density that determines the calculations of the formfactors in Eq. (2.21), i.e. the form factors all have even parity, see Fig. 15. Therefore the frequency shifts will carry no information about the odd part of the temperature perturbation, There are three approximate reflection symmetries in the LIGO test masses that result in degeneracy. Symmetry front to back z → −z, the symmetry right to left x → −x and symmetry up and down y → −y. The degeneracy associated with these symmetries result in two thermal profiles that are related by one of these symmetries, being indistinguishable. As a practical illustration, in Figure 4 we show two different thermal profiles as well as the change in mode frequencies computed in COMSOL for each of the profiles. The thermal profiles have intentionally been selected to have symmetry right to left. They are both 2D gaussian profiles across the mirror surface, uniform in depth. As expected the changes are almost equal, with precision of approximately 1%.
With the approximate symmetries of the Advanced LIGO test mass the maximum information that can be inferred from the frequency shifts is the symmetrized temperature distribution defined on one octant of the test mass:  for x, y, z > 0, here Σ denotes the summation over all possible combination of signs of x, y, z and the origin is assumed to be located at the center of mass of the test mass.
There may be a way of breaking some of the degeneracy by measuring other temperature-sensitive observables such as the distortion of the mirror's surface, or the thermal lensing of light passing through the test mass, for example with the Hartmann sensor. However, we do not consider these possibilities any further and leave their consideration for future work.

C. 3-dimensional temperature reconstruction using Singular Value Decomposition
Suppose now that we are considering properly symmetrized temperature fields so that f i (r) do form a complete basis. We should still exercise caution using Eq. (4.4) for the temperature field reconstruction. Similarity between some formfactors f i (r) means the C matrix is ill-conditioned (one or more eigenvalues are close to zero), the inversion becomes numerically unstable, leading to unreliable results. A common way of dealing with ill-conditioned matrices is to regularize the matrix by singular-value decomposition (SVD). The conversion matrix is decomposed into orthogonal matrix U, diagonal matrix C ′ and another orthogonal matrix V.
We adopt the convention that C ′ is defined with values sorted from largest to smallest along the diagonal.
Since C is real, the Hermitian transpose can be replaced by a regular transpose V * = V T . By removing eigenmodes associated with small eigenvalues, we reduce the dimensionality of C ′ to N − α by removing the α smallest elements of the complete diagonal matrix C ′ and the α associate eigenvectors in U and V. If α is suitably chosen, the resulting "regularized" matrix is numerically invertible. In what follows an example of singular value decomposition applied to simulated eigen-frequency data is demonstrated. The form factors for the first 225 eigenmodes of the test mass are calculated in COMSOL, each form factor is defined by the 60000 vertex elements existing in the three dimensional domain of the test mass. The inner product defined in Eq. (2.2) is performed to determine the conversion matrix C. Then the singular value decomposition is performed with Eq. (4.9). The relative numerical value of the eigenvalues (diagonal elements of C ′ ) provides a measure of the additional information that can be recovered by adding each new element in the SVD. These values are plotted in Figure 5. From the figure it can be seen that using more than 100 SVD elements does not provide a significant increase in information. It is also interesting to consider the shape of the largest elements in C ′ as these represent the temperature distribution components that will be most easily recovered. A selection of C ′ eigenfunctions are shown in Appendix B.
To demonstrate the usefulness of SVD, we consider a rotationally symmetric temperature distribution where r is the radial cylindrical coordinate and T 0 is a constant, and compute using Eqs. (2.1) and (2.21) the , without any significant difference in quality. However if we now assume that the eigenfrequency measurements are not perfect and contain errors, we note a marked difference in the quality of reconstructed temperature fields. As an example, we add a random frequency error drawn from a normal distribution with width 0.1 mHz to each analytically computed eigen-frequency change, and then compute the temperature fields from this erroneous data set. We observe that the truncated SVD inversion in Figure  6  the inversion that uses matrix C directly in Figure 6 (c). The latter is distorted due to errors in poorly resolved eignmodes. The optimal choice for α, the number of excluded eigenvectors, can be estimated for any particular temperature distribution, with a particular noise distribution by comparing the rms error of the temperature estimate b rms = V |dT − δT |dr. (4.11) In figure 7, the example plot of rms error as a function of the number of SVD elements used is shown. This example uses the same data as Figure 6. The optimal number of SVD elements is 86. The reconstructed temperature field with 86 elements is shown in Figure 6 d.

V. REALISTIC TEMPERATURE DISTRIBUTIONS
It is useful to reconstruct a temperature field that does not possesses all the symmetries of the test mass in order to see how our analysis works in real-world conditions. In this section we demonstrate that the recovery of symmetrized temperature distribution is indeed possible, circumventing the completeness problem due to test mass symmetries. If the temperature distribution is not symmetric in the same manner as the test mass symmetry the resulting rms error of the temperature distribution is large. In Fig. 9 this is demonstrated. In this case a a 100 kW beam on an optic with a uniform 1 ppm coating absorption is simulated resulting in a temperature distribution that is relatively higher on the high reflectivity surface and relatively cooler on the opposing surface ( Figure 8 top left and right panels respectively).
The estimated temperature distribution from the change in eigen-frequencies has roughly the same temperature distribution on the high reflectivity surface and This can be appreciated by comparing the rms error of the total test mass temperature distribution (blue) and the rms error of a model that uses half the test mass averaged with a reflection symmetry in the z-axis in (red) Figure 9, i.e., δT sym (x, y, z) = (1/2)ΣδT [x, y, ±z], where z > 0. More generally the average temperature distribution over one octant of the test mass may be computed when considering an arbitrary temperature distribution. While some information is lost, the symetrized temperature distribution still provides useful information. One potentially useful example is the measurement of the radial position of a beam on a test mass.

VI. PARAMETER ESTIMATION USING THE 3D TEMPERATURE FIELD
In the previous sections it was demonstrated measurements of a set of eigen-frequencies can be used to measure temperature distribution. The temperature distribution in the optic at LIGO may be defined by a relatively small number of parameters [12] of a thermal model. This limited model is described in Table II. In this section we show that the measurements of eigenfrequecnies can be used to measure specific thermal model parameters. We demonstrate how this can be done with Advanced LIGO data using a Bayesian approach. Eigenfrequency information can be collected during normal Advanced LIGO operation. We note that thermal conductivity affects the time evolution of the temperature inside the test mass. Therefore rather than using eigenfrequency measurements at one point in time we should use measurements over a time span δω i (t j ).
Some parameters such as laser power absorbed in the mirror coating (from the previous section) affect the temperature distribution in a linear fashion: T (P abs ) ∝ P abs × T (P abs = 1). (6.1) In this case, the power absorbed in the coating P abs may be inferred by linear regression: Previously, a single eigenmode has been used in such an analysis [12]. Linear regression using many eigenmodes benefits from more data and therefore less susceptibility to noise. Using more than one eigenmode also provides additional robustness against errors in different thermal model parameters. As errors produce temperature distributions with components orthogonal to the temperature distribution of interest, only the temperature distribution component common to both model parameters will affect the result. For a concrete example, consider a numerical experiment with P abs = 0.2 ppm. Consider that there is a 10% error in thermal conductivity such that k = 1.38 for the calculation of δω i (t j ) and k = 1.52 for the calculation of δω i (t j |P abs ). We then compareP abs calculated with Equation 6.2 with one eigenfrequency (the 6 th mode at 9330 Hz) and 100 eigenfreuencies (ranging from 5740 to 24888Hz). With one eigenfrequency, the estimate is biased and inaccurateP abs = -0.162±1.5 ppm. Using 100 eigenfrequencies, the estimate is precise and accuratê P abs = 0.199±0.005 ppm. This demonstrates the significant improvement in accuracy and robustness achieved by using a large number of eigenfrequencies.
More generally a set of thermal model parameters Γ may be estimated by locating the peak in the likelihood function.
3) In this paper we explore the likelihood function over various parameter spaces to determine what information is most easily recovered using this technique.
Transients in temperature are caused by laser light being absorbed in the test-mass mirror coating, changes in ring heater power and changes in ambient temperature. In this section, we focus on transients caused by laser light absorbed in the mirror coating as this is the most common thermal transient in LIGO optics. The thermal model for such a transient in its simplest form is defined by the optic geometry, the material properties of fused silica, and the properties of heat sources defined in Tables I and II. Laser light is absorbed in the mirror surface. Thermal equilibrium is attained when the radiative cooling to the thermal bath balances the heat load on the mirror surface.
Information recovered from multiple eigen-frequencies represent temperature gradients in the optic. Thermal gradients dissipate on a time scale proportional to the thermal gradient length scale. Therefore, the timescale of interest depends on characteristic length scale of the expected temperature field. For illustrative examples, and to keep computational costs low, we assume the eigenfrequneices are measured 10 times, with t i logarithmically spaced between 3 seconds and 10 hours.
A simulated example transient of a selection of eigenfrequencies is shown in Figure 10 for two different values of thermal conductivity k; the difference between the eigenfrequencies evaluated with different thermal conduc- tivity is shown as a dot-dash line. Note that most of the action, where mode frequencies change relative to each other, happens between a few hundred and a few thousand seconds. This is therefore the region we would expect to get most information regarding differences in temperature distribution for different thermal conductivity.
The log likelihood function is calculated for a simulated measurement point of k = 1.381 and is plotted in Figure 11. The log likelihood function peaks at the simulated measurement point, showing that this method can be used to infer properties like the thermal conductivity. This example assumes no uncertainty on any other parameters.
As apparent from Table II there are many model parameters that are subject to significant uncertainties. To get a sense of what parameters may be constrained using the method defined in this paper we investigated parame-ters in pairs. In Figure 12, the likelihood function is plotted over the parameter space of absorbed power P Abs and thermal conductivity k. The injected point is marked red. It can be seen that with small additional noise (1 nHz) other than quantization noise of the finite element simulation, both parameters are well constrained (coloured contour lines). However with 0.1 mHz measurement noise, the absorbed power is well constrained while the thermal conductivity can not be well constrained (grey lines with values indicated).
FIG. 12: likelihood function of point absorber power and thermal conductivity with minimal measurement noise (1 nHz) (coloured contour lines) and with 0.1 mHz measurement noise (gray contour lines). It can be seen that with measurement noise the absorbed power may be constrained while the thermal conductivity may only be minimally constrained These simulation were done for many pairs of parameters in Table II. Generally the absorbed power and the beam radial position are reasonably well constrained while the X and Y position estimates are less well constrained. Emissivity can be constrained in a similar manner to thermal conductivity. Other parameters are less well constrained.
Finally, in this section we show how this technique can be used to estimate thermal model parameters that are not accessible with Hartmann wavefront sensors. The thermal model of Table II assumes uniform absorption in the mirror high reflectivity coating. This model has recently been shown to be inadequate [3]. Point absorbers on the high reflectively surface of the test mass produce significant heating. The position of such a point absorber can be recovered well with the methods presented here. However this information is also accessible with the Hartmann wavefront sensor. In the following simulation we imagine a situation where instead of coating point absorbers there is a point absorber in the bulk of the test mass. The point absorber is a 30 um, 10% absorption feature. Figure 13 shows the likelihood function for the data given the point absorber location along with the true value of the point absorber location in red. While the distribution is bimodal, in this particular case, the absorber position is recovered as the maximum in the likelihood function, however with realistic measurement noise a bias is introduced. The thermal transient due to change in laser power is a common occurrence happening about once per day. Therefore a multi-parameter estimation may be arbitrarily refined using a Bayesian approach where the posterior distribution of the thermal model parameters inferred from one transient in laser power is used as the prior distribution for the subsequent measurement.

VII. CONCLUSION
Establishing robust thermal control of the test masses is one of the important tasks that will allow LIGO and Virgo to attain design goals. In this paper we provided an efficient method of computation of the vibrational mode frequency response to a temperature perturbation in the test mass. We demonstrated that the method may be inverted, enabling the conversion of vibrational mode frequency measurements into temperature distribution information. Finally, it was demonstrated that parameters of the test-mass thermal model may be estimated with improved precision using this temperature distribution information. Symmetries of the test mass prevent the recovery of complete 3D temperature distribution information, only symmetric components of the temperature distribution may be recovered. In principle, information from the Hartmann sensor could be used to break degeneracy between these symmetries and provide more information on the 3D temperature distribution. The framework described in this paper is demonstrated to provide useful coating absorption estimates and may allow estimates of several other thermal model parameters. However, this is dependent the nature of the measurement noise. Further experimental work and on-site measurements are needed to determine how the techniques proposed in this paper will be helpful for thermal control of the test masses.
The initial stages of this research were supported by Levin's Australian Research Council (ARC) Future Fellowship and Blair's PhD scholarship. Later work was supported by Dr. Blair's Caltech postoctoral fellowship and ARC DECRA DE190100437. ET is supported through ARC Future Fellowship FT150100281 and ARC Centre of Excellence CE170100004. Figure 14 shows a pair of eigenfrequencies that have very similar form factors. These similarities make the conversion matrix rank deficient and thus singular value decomposition is required. Figure 15 shows a selection of vibrational wavefunctions and their associated form factors. Blue regions of the form factors indicates areas of the test mass where the particular mode frequency is insensitive to temperature variation. Red regions are areas where the mode frequency is sensitive to temperature variation. These eigenfrequencies are measured in LIGO data and therefore formfactors represent temperature distribution spatial scale factors that should be measurable.