High-accuracy numerical models of Brownian thermal noise in thin mirror coatings

Brownian coating thermal noise in detector test masses is limiting the sensitivity of current gravitational-wave detectors on Earth. Therefore, accurate numerical models can inform the ongoing effort to minimize Brownian coating thermal noise in current and future gravitational-wave detectors. Such numerical models typically require significant computational resources and time, and often involve closed-source commercial codes. In contrast, open-source codes give complete visibility and control of the simulated physics, enable direct assessment of the numerical accuracy, and support the reproducibility of results. In this article, we use the open-source SpECTRE numerical relativity code and adopt a novel discontinuous Galerkin numerical method to model Brownian coating thermal noise. We demonstrate that SpECTRE achieves significantly higher accuracy than a previous approach at a fraction of the computational cost. Furthermore, we numerically model Brownian coating thermal noise in multiple sub-wavelength crystalline coating layers for the first time. Our new numerical method has the potential to enable fast exploration of realistic mirror configurations, and hence to guide the search for optimal mirror geometries, beam shapes and coating materials for gravitational-wave detectors.


Introduction
Brownian coating thermal noise is the limiting noise source for current-generation, ground-based gravitational-wave detectors in their most sensitive frequency bands. For instance, following the A+ upgrade anticipated for completion in the mid 2020s, the Laser Intererometer Gravitational-Wave Observatory (LIGO) detector noise is dominated by Brownian coating thermal noise at frequencies f ∼ 100 Hz [1]. This noise arises from thermal fluctuations in the reflective coatings of the detectors' test masses [2].
Therefore, a reduction of the Brownian coating thermal noise directly increases a detector's sensitivity and thus its astronomical reach. Theoretical models of Brownian coating thermal noise are important for working toward this goal. Thermal noise modeling typically follow the approach pioneered by Levin [3], which computes the thermal noise in terms of an auxiliary elasticity calculation using the fluctuationdissipation theorem [4,5,6]. While an approximate analytic solution is well known in the limit where coating thickness and edge effects can be neglected, numerical calculations of thermal noise are necessary to study effects that arise from the finite test-mass size, the finite coating thickness, and from crystalline materials.
In this article we calculate Brownian coating thermal noise by numerically solving the auxiliary linear elasticity problem. Such numerical simulations typically adopt a conventional finite-element approach, as some of the authors did in Ref. [7]. These methods are widely used, but achieving high accuracy with them can require significant computational resources and time, because of their relatively slow rates of convergence.
For the first time to our knowledge, we apply a discontinuous Galerkin (DG) method to model Brownian coating thermal noise. DG methods are well suited to this problem because they can retain high-order convergence in the presence of discontinuities, which arise at the interfaces between the mirror substrate and its reflective coatings. In this article, we extend the DG method for elliptic equations presented in Ref. [8] to problems with discontinuous material properties. With this extension, our method converges exponentially with resolution, allowing us to solve coating thermal noise problems numerically at high accuracy using considerably less computational resources and time than conventional finite-element methods.
We implement the numerical method and the elastostatic equations in SpECTRE [9], a new open-source numerical relativity code. While SpECTRE's primary aim is to model merging black holes and neutron stars, the elliptic solver needed to construct initial data for such simulations is also very well positioned to solve the DG-discretized elastostatics equations for thermal noise modeling [10,8]. As an open-source code, our approach has advantages compared to the closed-source and commercial solutions that are often adopted: we can directly control the physics incorporated in the calculation, we can assess the accuracy and convergence rate of our simulations in a straightforward way, and our results are reproducible with publicly available software. Our code also benefits from SpECTRE's task-based parallelism approach, implemented using the Charm++ [11] library, enabling our code to efficiently scale to large numbers of compute cores [12].
This article is organized as follows. Section 2 summarizes the elastic problem to be solved and presents the discontinuous Galerkin numerical method. Section 3 presents our results using this method to model thermal noise in cylindrical mirrors with thin coatings. We discuss our results and future work in Sec. 4.

Methods
In this section, we formulate the auxiliary elasticity problem based on Refs. [3,7], discretize it with the discontinuous Galerkin scheme developed in Ref. [8], and outline the numerical method we employ to solve the discretized problem with the SpECTRE code [10]. Section 2.2.3 details a novel extension of this method to handle discontinuous material properties at layer interfaces.

Auxiliary elasticity problem
We consider a gravitational-wave detector that measures the position of a test mass with a laser beam with a Gaussian intensity profile Here, r is the cylindrical radial coordinate from the center of the beam with width r 0 . The intensity profile is normalized so that The laser beam effectively measures a weighted average q of the displacement Z of the test mass surface, As shown by Levin [3], Brownian thermal noise can be calculated from the energy dissipated in an auxiliary elastic problem. Specifically, to compute the thermal noise at frequency f , one applies an oscillating pressure to the face of the mirror with frequency f , with a pressure distribution profile p(r) equal to the beam intensity, and with an amplitude F 0 . In this auxiliary problem, the energy W diss will be dissipated in each cycle of the oscillation. The fluctuation-dissipation theorem relates this dissipated energy W diss to the thermal noise, specifically to the power spectral density S q associated with q, 1 where T is the mirror temperature and k B is Boltzmann's constant. Because W diss ∝ F 2 0 , it follows that S q does not depend on the overall amplitude F 0 . For frequencies f ∼ 100 Hz much lower than the resonant frequencies f ∼ 10 4 Hz of the test-mass materials, the dissipated power can be computed using the quasistatic approximation. In this approximation, a static pressure is applied to the mirror with amplitude F 0 and profile p(r), and the dissipated energy can be written as where U is the potential energy stored in the deformation of the test-mass and φ is the material's loss angle determined by the material's imaginary, dissipative elastic moduli.
Therefore, our goal in this article is to solve the equations of elastostatics for the deformation of the test mass, when its surface is subjected to an applied pressure with profile p(r). Here, T ij is the stress and we adopt the Einstein summation convention so that repeated tensor indices are summed over. The source f j is the force density acting on each volume element of the mirror as a function of position x, which vanishes in our situation, f j = 0. The pressure acting on the external surface of the test mass will be represented by suitable boundary conditions. Equation (6) is an equation for the displacement vector field u i (x), which describes the deformation of the elastic material as a function of the undeformed coordinates. The symmetric part of the gradient of the displacement vector field is the strain For sufficiently small F 0 , the strain is proportional to the applied stress, where the constitutive relation Y ijkl (x) captures the elastic properties of the material in the linear regime. The constitutive relation is symmetric on its first two indices, on its last two indices, and under exchange of the first pair of indices with the second pair of indices. Inserting Eqs. (7) and (8) into Eq. (6) yields the equations of linear elasticity, which we will solve numerically. We consider materials with either amorphous or cubic-crystalline constitutive relations. The coating may consist of only a single material, or of multiple layers with different materials. The amorphous constitutive relation is isotropic and homogeneous, with Lamé parameter λ and shear modulus µ. 2 A cubic-crystalline material is characterized by the constitutive relation where c 11 , c 12 and c 44 are three independent material parameters. The constitutive relation in Eq. (9) is composed by discontinuously choosing either Eq. (10) or Eq. (11) in each layer of the material. After solving the linear elasticity equations, Eq. (9), the potential energy is evaluated by an integral over the volume of the material, For a material with a thin, reflective coating with different elastic properties than the substrate, the dissipated energy, Eq. (5), decomposes as [14] where U sub and φ sub are the potential energy and loss angle of the substrate, respectively, while U coat and φ coat are the potential energy and the loss angle of the coating. Note that a material can also have different loss angles associated with the different independent elastic moduli of a material [15]. We do not consider further decompositions of the elastic potential energy in this article, but note that such quantities can straightforwardly be extracted from our simulations. The different loss angles only affect the computation of the thermal noise by Eq. (4) from the elastic potential energy extracted from our simulations, but they do not affect our simulations or our numerical method in any way. An approximate analytic solution exists for amorphous materials in the limit where the coating thickness d is small compared to both the size of the mirror and the width r 0 of the pressure profile. The approximate coating thermal noise is 3

Discontinuous Galerkin discretization
We employ the discontinuous Galerkin (DG) scheme detailed in Ref. [8] to discretize the elasticity problem, Eq. (9). We summarize the discretization scheme in this section, and extend it to problems with discontinuous material properties.

Domain decomposition
We simulate a cylindrical mirror in three dimensions with radius R and height H. The cylinder axis coincides with the z axis of our coordinates, and the plane z = 0 represents the surface of the mirror on which the external pressure p(r) is applied. We decompose the cylindrical domain Ω = [0, R] × [0, 2π) × [0, H] into a set of nonoverlapping elements Ω k ⊂ Ω shaped like deformed cubes, as illustrated in Fig. 1a (h refinement). Each element carries a coordinate map from the Cartesian coordinates x ∈ Ω k , in which the elasticity equations (9) are formulated, to logical coordinates ξ ∈ [−1, 1] 3 representing the reference cube, as illustrated in Fig. 1b. The coordinate map to the reference cube is characterized by its Jacobian, with determinant J and inverse (J −1 ) On the reference cube we choose a set of N k,i Legendre-Gauss-Lobatto (LGL) collocation points in each dimension i (p refinement).
Fields are represented numerically by their values at the collocation points. We denote the set of discrete values for the displacement vector field u i within an Top: Geometry of our layered cylindrical domain, with the laser beam indicated in red. Four wedge-shaped elements envelop a cuboid. Another set of wedges extends to the outer radius of the cylinder. In z direction the cylinder is partitioned into layers that can have different material properties (black and gray). The substrate layer has a logarithmic coordinate map in z direction and is split in two twice in this example (thin horizontal lines). Bottom: The coordinate transformation ξ(x) maps an element to a reference cube ξ ∈ [−1, 1] 3 with logical coordinate-axes ξ = (ξ, η, ζ). In this example we chose N k,ξ = 3 and N k,η = 4 LGL collocation points along ξ and η, respectively.
and the collection of discrete displacement vector field values over all elements as u i . The values at the collocation points within an element define a three-dimensional Lagrange interpolation, where the basis functions ψ p (ξ) are products of Lagrange polynomials, based on the collocation points in the three logical directions of the element Ω k . Since Eqs. (17) and (18) are local to each element, fields over the entire domain are discontinuous across element boundaries.

DG residuals
To formulate the elasticity equations in first-order form for the DG discretization, we use the symmetric strain S kl as auxiliary variable. Following Ref. [8], we first compute the discrete auxiliary variables on the computational grid as where we make use of the discrete differentiation matrix D i := M −1 MD i , the mass matrix the stiffness matrix the lifting operator and L := M −1 ML on the element Ω k [8]. The integral in Eq. (22) is over the boundary of the element, ∂Ω k , where n i is the outward-pointing unit normal one-form and J Σ is the surface Jacobian. The symbol · emphasizes matrix multiplication with the field values over the computational grid of the element. In a second step, we compute the DG residuals in strong form [8], which represent the set of algebraic equations for the values u i of the displacement vector field on the computational grid that we solve numerically.

Numerical flux
The quantities (n (k u l) ) * and (n i Y ijkl S kl ) * in Eq. (23) denote a numerical flux that couples grid points across nearest-neighbor element boundaries. We employ the generalized internal-penalty numerical flux developed in Ref. [8], with one notable extension. Contrary to Ref. [8] we allow neighboring elements to define different constitutive relations, meaning Y ijkl (x) can be double-valued on shared element boundaries. 4 Therefore, we define the quantity where "int" denotes the interior side of an element's shared boundary with a neighbor, and "ext" denotes the exterior side, i.e. the neighbor's side. With this quantity we can define the numerical flux where n ext i = −n int i for the purpose of this article. Equation (25) is the generalized internal-penalty numerical flux defined in Ref. [8], with a choice between Y ijkl int , Y ijkl ext , and Y ijkl * for every occurrence of the constitutive relation. The particular choice in Eq. (25) ensures that the numerical flux remains consistent, meaning that i T ij and n int (k u int l) = −n ext (k u ext l) . In particular, note that the penalty term in Eq. (25b) vanishes when the displacement is continuous across the boundary, and that the numerical flux admits solutions where the stress is continuous across the boundary but the strain is not. Such solutions may arise in a layered material under stress, because the layers remain "glued together" but each layer responds to the stress differently.
The penalty function in Eq. (25b) is where we make use of the polynomial degree p and a measure of the element size, h, orthogonal to the element boundary on either side of the interface, as detailed in Ref. [8]. We choose C = 100 in this article.

Boundary conditions
We impose boundary conditions through fluxes, i.e. by a choice of exterior quantities in the numerical flux (25). Specifically, on external boundaries we set where we choose either u i b to impose Dirichlet boundary conditions, or n int i T ij b to impose Neumann boundary conditions on the boundary collocation points, and set the respective other quantity to its interior value.
For the thermal noise problem we impose the pressure induced by the laser beam, as Neumann boundary condition on the z = 0 side of the cylindrical mirror, where p(r) is the laser beam profile given in Eq. (1). On the side of the mirror facing away from the laser we impose as Dirichlet boundary condition, and on the mantle we impose as Neumann boundary condition. Equation (29) means that the back of the mirror is held in place, whereas Eq. (30) implies no pressure on the sides, which however, are free to deform in response to the pressure applied to the front.

SpECTRE elliptic solver
Once discretized, the linear algebraic equations (23) are solved numerically for the displacement vector field values u i on all elements and grid points in the computational domain. As is typical for discretized elliptic equations, Eq. (23) defines a matrix equation where u denotes the set of all N DOF = 3 × N points = 3 × k N k displacement vector field values in the computational domain, and A is a matrix with N DOF ×N DOF entries. To solve Eq. (31) means inverting the matrix A. However, as the resolution of the computational domain increases, the matrix A easily becomes too large to construct explicitly, to store on an ordinary computer, and to invert directly. Therefore, we solve Eq. (31) with the elliptic solver component of the open-source SpECTRE code [9,10]. It employs an iterative generalized minimal residual (GMRES) algorithm to solve Eq. (31) to the requested precision. A multigrid preconditioner accelerates the GMRES algorithm by supporting each iteration with an approximate solution from a hierarchy of successively coarser grids. On every grid, an additive Schwarz smoother decomposes the problem into many overlapping subproblems, one per element in the domain, which are solved independently and in parallel. The subproblems are distributed across the cores of a computing cluster by a task-based parallelization paradigm. The elliptic solver is described in detail in Ref. [10].

Results
Our simulations with SpECTRE were performed on one or more 16-core compute nodes, each with 64 GB of memory and two eight-core Intel Haswell E5-2630v3 processors clocked at 2.40 GHz, connected with an Intel Omni-Path network. We distribute the elements that compose the computational domain evenly among cores, leaving one core per node free to perform communications.
We compare our results to previous work using an open-source finite element code to calculate the Brownian coating thermal noise for amorphous and crystalline materials. Its methods are described in Secs. 2.4-2.6 of Ref. [7]. The code was built using the deal.ii [16,17] finite element framework and we henceforth refer to it as deal.ii. It adopted a standard weak form of the elastostatic equations, discretized them using a conventional finite element approach, and solved them using deal.ii with the PETSc [18] conjugate gradient linear solver and the ParaSAILS preconditioner in the hypre [19] package. The deal.ii code relies on the Message Passing Interface (MPI) for parallelization.

Single-coating comparison
First, we consider the single-coating scenario investigated in Ref. [7] and demonstrate the superior performance of our new approach. We choose the parameters listed in Ref. [7], Table 1 for a cylindrical mirror of radius R = 12.5 mm with a single d = 6.83 µm thin effective-isotropic AlGaAs coating. We simulate the scenario both with the deal.ii approach employed in Ref. [7] and with our new approach with the SpECTRE code. Figure 2 presents the numerical precision and computational cost of both approaches. To assess the numerical precision we successively increase the resolution in both codes. In SpECTRE we increase the resolution by incrementing the number of grid points in all dimensions of all elements in the domain by one, and in deal.ii we employ an adaptive mesh-refinement scheme [7]. We compute the error in the elastic potential energy relative to a high-resolution reference configuration simulated in SpECTRE, for which we have split all elements of the highest-resolution configuration included in Fig. 2 in half along all three dimensions (see also Fig. 1a). Hence, the highresolution reference configuration has ∼ 79 grid points per dimension.
We find that both codes converge to the same solution, but our new approach in SpECTRE achieves about four orders of magnitude higher accuracy than the deal.ii approach using the same number of grid points. Furthermore, our new approach simulates this scenario with sub-percent error in only 30 s on 15 cores, for which the deal.ii approach required multiple hours on 324 cores. Our new approach also achieves a fractional error below 10 −5 in only half a core-hour, or two minutes of real time, which was prohibitively expensive with the deal.ii approach.

Accuracy of the approximate analytic solution
Second, we study the accuracy of the approximate analytic solution for the singlecoating thermal noise, Eq. (14), using the superior numerical precision we can now achieve over the results presented in Ref. [7]. The approximate solution holds for a thin coating, d/r 0 1, a semi-infinite mirror, r 0 /R 1 and d/R 1, and for isotropic-homogeneous materials. Therefore, it does not capture the finite-size effects included in our simulations, and approximates the crystalline AlGaAs coating as an amorphous material.
To assess the magnitude of the finite-size effects, we employ the simulations detailed in Sec. 3.1, which use the same effective-isotropic model for the AlGaAs  Thermal noise in an AlGaAs-coated mirror computed from the approximate analytic solution (14) and from our numerical simulations. The effective-isotropic simulation (black) retains the amorphous approximation for the material, but includes finite-size effects. The crystalline simulation (red) eliminates this approximation. Previous simulations with the deal.ii approach are shown in lighter colors to the left.
coating that underpins the approximate analytic solution. Figure 3 presents both the thermal noise computed from the simulations and the approximate analytic solution (black). Error bars are computed as ∆ S coat q / S coat q = 1/2 ∆U coat /U coat from the relative numerical error in the elastic potential energy. While Ref. [7] estimated the magnitude of finite-size effects for this problem to 7 %, we can now report that their simulations captured the effect to 7.5(2) %. With our new numerical method, we can make this statement more precise and report a finite-size effect of 7.616 649(6) %.
To assess the magnitude of the amorphous approximation to the crystalline coating material, we repeat the simulations with a crystalline constitutive relation and the parameters listed in Ref. [7], Table 1. The thermal noise computed from these simulations is presented in Fig. 3 as well (red). We refine the estimate of 4 % from Ref. [7] to 4.5(2) %, and report 4.667 990(6) % using our new numerical method.
Note that we report only numerical errors from our simulations here. However, the parameters that define our simulations, such as coating and substrate material properties, are typically measured experimentally and carry significant uncertainties. For example, elastic moduli recently reported in Ref. [20] were measured on the percent level. Computing the thermal noise by Eq. (4) also involves loss angles φ which can be measured on the percent level as well [20]. Therefore, computational resources are spent most effectively to drive numerical errors below sub-percent levels and no further. With our new numerical methods we achieve sub-percent accuracy with a fraction of the computational resources required before (see Sec. 3.1).

Multiple sub-wavelength crystalline coatings
Finally, we apply our new computational approach to a scenario that presents many of the challenges we expect for applications to realistic mirror configurations. We  Elastic variables along the z-axis for a mirror with multiple thin coating layers, simulated with SpECTRE. Gray layers are crystalline AlGaAs, and white regions are fused silica. The nine coating layers have a combined thickness of 6.83 µm and the material extends to z = 12.5 mm outside the range of this plot. simulate a cylindrical mirror of the same radius R = 12.5 mm as before, but split the d = 6.83 µm thin coating into nine layers, so the thickness of each coating layer is below the typical 1 µm wavelength of the laser. The coating layers alternate between fused silica and crystalline AlGaAs, with the elastic moduli c 11 , c 12 and c 44 listed in Ref. [7], Table 1. Neither sub-wavelength coatings nor multiple layers were simulated in Ref. [7], but our new computational approach in SpECTRE achieves both. Figure 4 presents our numerical solution of this scenario. Our new computational approach based on discontinuous Galerkin methods resolves the thin coating layers at high accuracy without spurious oscillations. Figure 5 presents the numerical precision of the solution. We increase the resolution by incrementing the number of grid points per element and dimension and compute the relative error to a high-resolution reference configuration, as we did in Sec. 3.1. The error converges exponentially, which is a feature of our discontinuous Galerkin method with grid boundaries placed at the layer interfaces.

Discussion
We have presented a new numerical method to model Brownian thermal noise in thin mirror coatings based on a discontinuous Galerkin (DG) discretization. With our new method, we model thermal noise in a one-inch cylindrical mirror with a microns-thick coating at unprecedented accuracy at a fraction of the time needed in a previous, conventional finite-element approach [7]. Using these high-accuracy simulations, we Relative error |U − Uref
Convergence test of the elastic potential energy in each coating layer, and in the substrate. Our new numerical method achieves exponential convergence despite the discontinuous material properties.
find that a commonly-used approximate analytic solution overestimates the coating thermal noise for this problem by 7.6 % when taking only finite-size effects into account, and by 4.7 % when modeling it as a crystalline material, which refines a previous estimate in Ref. [7]. We also demonstrate that, unlike the approach in Ref. [7], our new method is capable of resolving multiple sub-wavelength coatings, including coatings of a cubic-crystalline material. Our new numerical method is implemented in the opensource SpECTRE code and the results presented in this article are reproducible with the supplemental input-file configurations.
We found that it is crucial for the success of our new method that the interfaces between layers of different materials coincide with element boundaries in our computational domain. Then, our discontinuous Galerkin discretization with a suitable choice of numerical flux converges exponentially, achieving high accuracy with a small number of grid points. The scheme can potentially be improved in future work. Most notably, an adaptive mesh-refinement (AMR) algorithm would have great potential to further improve the accuracy and efficiency of the scheme, by distributing the resolution in the computational domain to regions and dimensions where it is most needed.
Furthermore, the elliptic solver in the SpECTRE code that we employ to solve the discretized problem numerically can be improved to accelerate thermal-noise calculations. The calculations we have presented in this article require a few hundred solver iterations to converge, or up to ∼ 1400 for our highest-resolution simulation with multiple sub-wavelength crystalline coatings. While simple configurations complete in seconds or minutes of real-time on 15 cores, where the previous approach needed hours on 324 cores, the more challenging configurations, which were prohibitively expensive with the previous approach, solve in about an hour on 45 cores.
We expect additional speedup with further improvements to the elliptic solver algorithm in SpECTRE. In particular, improvements to its multigrid preconditioner have great potential to speed up the simulations. The multigrid algorithm relies on solving the problem approximately on coarser grids to resolve large-scale modes in the solution. It currently cannot coarsen the grid any further than the size of each coating layer because the layers define the material properties. To accelerate the calculations, we intend to let the multigrid algorithm combine layers with different materials into fiducial coarse layers with effective material properties. This approach is possible because the partitioning of the domain into layers is necessary only to define material properties, not to define the cylindrical shape of the domain. Reference [10] shows that the multigrid algorithm can achieve resolution-independent iteration counts when the domain can be coarsened sufficiently. Note that the fiducial coarse layers affect only the convergence speed of the solver and do not change the solution once the solver has converged.
Our numerical models of thermal noise have the potential to inform upgrades that increase the sensitivity of gravitational-wave detectors, using the advanced computational technology that we develop for numerical-relativity simulations in the SpECTRE code. In the future, we intend to apply our new numerical method to simulate Brownian thermal noise in more realistic mirror configurations and materials that are under consideration for current and future gravitational-wave detectors, such as the optimized configuration found in Ref. [21]. While approximate analytic solutions can provide useful estimates, only numerical models can precisely quantify the finitesize effects of changing the mirror geometry. In particular, finite-size effects are more important for real gravitational-wave detectors than for tabletop experiments measuring thermal noise. Tabletop experiments often use small beam sizes to enlarge the thermal noise and hence make it easier to measure, whereas gravitational-wave detectors prefer large beam sizes to minimize thermal noise. Therefore, we plan to employ our new numerical method to explore realistic mirror configurations, with the goal of finding configurations that minimize Brownian coating thermal noise.