# Closed vortex state in 3D mesoscopic superconducting films under an applied
transport current

Leonardo Rodrigues Cadorim  Departamento de Física, Faculdade de Ciências,
Universidade Estadual Paulista (UNESP), Caixa Postal 473, 17033-360, Bauru-SP,
Brazil  Departement Fysica, Universiteit Antwerpen, Groenenborgelaan 171,
B-2020 Antwerp, Belgium  Lucas Veneziani de Toledo  Departamento de Física,
Faculdade de Ciências, Universidade Estadual Paulista (UNESP), Caixa Postal
473, 17033-360, Bauru-SP, Brazil  Wilson Aires Ortiz  Departamento de Física,
Universidade Federal de São Carlos (UFSCar), Caixa Postal 676, 13565-905, São
Carlos-SP, Brazil  Jorge Berger  Department of Physics and Optical
Engineering, Braude College, Karmiel, 2161002, Israel  Edson Sardella
Corresponding author<EMAIL_ADDRESS>Departamento de Física,
Faculdade de Ciências, Universidade Estadual Paulista (UNESP), Caixa Postal
473, 17033-360, Bauru-SP, Brazil

###### Abstract

By using the full 3D generalized time dependent Ginzbug-Landau equation we
study a long superconducting film of finite width and thickness under an
applied transport current. We show that, for sufficiently large thickness, the
vortices and the antivortices become curved before they annihilate each other.
As they approach the center of the sample, their ends combine, producing a
single closed vortex. We also determine the critical values of the thickness
for which the closed vortex sets in for different values of the Ginzburg-Ladau
parameter. Finally, we propose a model of how to detect a closed vortex
experimentally.

††preprint: APS/123-QED

##  I Introduction

One of the most outstanding physical phenomena in condensed matter theory is
the flux quantization in the Shubnikov phase of superconducting materials. In
bulk type-II superconductors, above a certain critical value of applied field,
magnetic flux penetrates the sample in the form of flux quanta (in units of
$\Phi_{0}=hc/2e$), the so called magnetic vortices. Superconducting vortices
are very important due to their possible use in many applications, ranging
from single photon detectors to quantum information. Moreover, the knowledge
of how vortices behave under various circumstances is fundamental to the
optimization of many electronic devices, since the vortex motion causes energy
dissipation in a superconductor.

In the core of a superconducting vortex, the order parameter that describes
the superconducting state vanishes and its phase changes by a multiple of
$2\pi$ when circulated around a closed loop that encloses the core. As the
vortex line tends to align with the direction of the magnetic field, in the
presence of an external field applied perpendicularly to a superconducting
film, the core of the vortex is rigorously a straight line with the currents
flowing around it. This is the so called Abrikosov vortex, first described in
his outstanding work [1].

However, the circularly shaped magnetic self-field induced by a current
produces an Abrikosov vortex with a ringlike shape, which is usually called a
closed vortex. This new type of solution was first found by Kozlov and
Samokhvalov [2] through the solution of the London equation and was
extensively studied further in Refs. [3, 4, 5, 6, 7, 8, 9, 10, 11, 12]. In
these works, the existence, dynamics and even stability of such closed
vortices in the presence of inhomogeneities were studied for unbounded
superconductors or superconducting samples in the shape of a cylinder. In the
latter case, due to matching symmetry between the field lines produced by the
current and the cylinder geometry, the vortex penetrates the superconductor
already in the form of a closed ring. Recently, these closed vortices were
also shown to exist in a more complex geometry, such as a superconducting
torus. [13]

In different scenarios, the formation of closed vortices was also recently
studied, for example, in Josephson junctions [14] submitted to an external
current, where it was shown that Josephson vortices can also be found in the
form of closed vortex loops and a procedure to their experimental observation
was introduced. Due to short life time of the Abrikosov closed vortices, their
experimental detection is a difficult task and has not been accomplished so
far. Recently, by using the microscopic theory of superconductivity Fyhn and
Linder [15] proposed an experimental setup to the observation of such objects
based on STM measurements. One of the most important results of the present
work is to provide a new procedure to the detection of closed vortex loops,
which is based on the measurement of the magnetic field profile produced by
them. It is worth mentioning that, in practice, half-closed vortices have been
demonstrated in RF cavities. [16]

In the present work, we study the physics of the closed vortices in a
different setup. Namely, we investigate how a closed vortex emerges and gets
annihilated in a superconducting slab under the presence of a transport
current. Unlike previous works on the subject, the geometry of our
superconductor does not favor the formation of a closed vortex. Here, this
object is created solely by the inhomogeneous action of the applied current in
different parts of the flux line. The study of the vortex dynamics is of
fundamental importance to the understanding of the resistive state of current
driven superconductors. Here, we show that, for appropriate thickness of the
slab and Ginzburg-Landau parameter $\kappa$, a closed vortex forms in the
process of annihilation of a vortex-antivortex (v-av) pairs induced in the
sample by the current self-field. As we show, the lines corresponding to the
vortex and the antivortex combine and form a closed loop when the pair is
annihilating. After this combination, the loop takes the form of a quasi-
ellipse, with the aspect ratio gradually decreasing until the collapse of the
loop. By increasing the thickness of the film, the loop tends to a circle.

The outline of this work is as follows. In Sec. II we present our model and
the formalism we used in order to solve the generalized Ginzburg-Landau
equations [17, 18]. In Sec. III we present and discuss the results obtained in
our simulations. Finally, we present our concluding remarks in Sec. IV.

##  II  Theoretical model

In this work, we rely on the generalized time dependent Ginzburg-Landau
(GTDGL) equation which is more suitable to describe the resitive state of
dirty superconductors in the non-equilibrium state [17, 18]. In dimensionless
units this equation is given by

$\displaystyle\frac{u}{\sqrt{1+\gamma^{2}|\psi|^{2}}}\left[\frac{\partial}{\partial
t}+\frac{1}{2}\gamma^{2}\frac{\partial|\psi|^{2}}{\partial t}\right]\psi=$
$\displaystyle=\left(\mbox{\boldmath$\nabla$}-i\textbf{A}\right)^{2}\psi+\psi(1-|\psi|^{2}),$
(1)

coupled with Ampère’s law

$\Sigma\frac{\partial\textbf{A}}{\partial
t}=\textbf{J}_{s}-\kappa^{2}\mbox{\boldmath$\nabla$}\times\textbf{h},$ (2)

where

$\textbf{J}_{s}={\rm
Im}\left[\bar{\psi}{(\mbox{\boldmath$\nabla$}-i\textbf{A})}\psi\right]$ (3)

is the superconducting current density.

Here, the temperature is in units of the critical temperature $T_{c}$; the
order parameter $\psi$ is in units of
$\psi_{\infty}(T)=\sqrt{\alpha(T)/\beta}$, where $\alpha$ and $\beta$ are two
phenomenological constants; the distances are measured in units of the
coherence length $\xi(T)$; the vector potential A is in units of $\xi
H_{c2}(T)$, where $H_{c2}$ is the upper critical field; the local magnetic
field $\textbf{h}=\mbox{\boldmath$\nabla$}\times\textbf{A}$ is units of
$H_{c2}(T)$; time is in units of the Ginzburg-Landau characteristic time
$\tau_{GL}=\pi\hbar/8k_{B}Tu$; the material dependent parameter
$\gamma=2\tau_{E}\Delta_{0}/\hbar$, where $\tau_{E}$ is the inelastic
electron-collision time, and $\Delta_{0}$ is the gap in the Meissner state;
the constant $\Sigma=4\pi\sigma D/c^{2}\xi^{2}(T)$, where $D$ is the diffusion
coefficient and $\sigma$ is the normal state electrical conductivity;
$\kappa=\lambda(T)/\xi(T)$ is the Ginzburg-Landau parameter, where
$\lambda(T)$ is the London penetration depth; and finally, the constant $u$ is
equal to 5.79, which is derived from first principles [17].

The original GTDGL equations take into account the scalar electrical potential
$\varphi$. Since they are invariant under the following gauge transformations

$\displaystyle\psi^{\prime}$ $\displaystyle=$ $\displaystyle e^{-i\chi}\psi,$
$\displaystyle\textbf{A}^{\prime}$ $\displaystyle=$
$\displaystyle\textbf{A}-\mbox{\boldmath$\nabla$}\chi,$
$\displaystyle\varphi^{\prime}$ $\displaystyle=$
$\displaystyle\varphi+\frac{\partial\chi}{\partial t},$ (4)

where $\chi$ is an arbitrary scalar function therefore, we conveniently use
the Weyl gauge [19] in which the scalar potential is constant and equal to
zero.

In the present work, we consider an infinite superconducting film carrying a
transport current as sketched in Fig. 1. The applied transport current is
introduced as follows. The superconductor is in the presence of an applied
electric field which is sustained by a DC transport current density which
flows along the $x$ axis, $\textbf{J}_{a}=J_{a}\hat{\textbf{x}}$ . In the
normal state, the vector potential has only the $x$ component. Thus, we
separate the vector potential and the local magnetic field into two
contributions, one coming from the normal state, and another one due to the
diamagnetic nature of the superconductor that tends to cancel out the field
induced by the applied current inside the sample. In other words, in Eqs.1-3
we substitute:

A $\displaystyle=$ $\displaystyle\textbf{A}_{0}+\textbf{A}_{1},$ (5) h
$\displaystyle=$ $\displaystyle\textbf{h}_{0}+\textbf{h}_{1},$ (6)

where A0 and h0 satisfy the following equations:

$\displaystyle\kappa^{2}\mbox{\boldmath$\nabla$}\times\textbf{h}_{0}=J_{a}\hat{\textbf{x}},\;\;\;\kappa^{2}\nabla^{2}A_{0x}=-J_{a}\,.$
(7)

Figure 1: Schematic view of the system under consideration: an infinitely long
superconducting sample of width $l_{y}$ and thickness $l_{z}$; only one unit
cell of length $l_{x}$ is shown. The transport current is applied in the $x$
direction. The encircling lines illustrate the line fields of the self-field
produced by the current. Two defects are introduced at the border of the
sample (black spots), in order to facilitate nucleation of v-av pairs.

The analytical solutions of equations (7) are given in the Supplementary
Material of Ref. [20].

Here, we solve the full 3D GTDGL equations numerically for an infinite
superconducting film of finite width and thickness (see Fig. 1). The
infinitely long film is divided into unit cells of dimensions
$(l_{x},l_{y},l_{z})$. We take into account the demagnetization effects.
Therefore, for numerical purposes, we must consider the unit cell inside a
simulation box (not shown in Fig. 1) of dimensions $(l_{x},L_{y},L_{z})$,
where $(L_{y},L_{z})$ are sufficiently larger than $(l_{y},l_{z})$ so that the
demagnetizing field $\textbf{h}_{1}$ vanishes far away from the
superconducting surfaces (for more details see Ref. [21]). In addition, the
normal components of the superconducting current density of must be zero on
the superconductor-vacuum interfaces. Then, the following boundary conditions
must be fulfilled:

$\displaystyle\hat{\textbf{n}}\cdot(\mbox{\boldmath$\nabla$}-i\textbf{A}_{0}-i\textbf{A}_{1})\psi$
$\displaystyle=$ $\displaystyle 0,\;\;\;\textrm{in}\;\partial\Omega_{sc},$ (8)
$\displaystyle\mbox{\boldmath$\nabla$}\times\textbf{A}_{1}$ $\displaystyle=$
$\displaystyle 0,\;\;\;\textrm{in}\;\partial\Omega,$ (9)

where $\partial\Omega_{sc}$ and $\partial\Omega$ stand for superconducting and
simulation box surfaces, respectively.

The two black spots in Fig. 1 represent two defects on the border of the
sample. They are introduced as an artifact to create an inhomogeneity in the
current, that induces the vortices and antivortices in the opposite sides of
the sample.

##  III Results and discussion

###  III.1 Parameters and Methodology

The above equations are discretized by using the standard link-variable method
as described in Ref. [22]. This algorithm is then implemented in Fortran 90
programming language and run in a GPU (Graphics Processing Unit) accelerated
forward-time-central-space scheme.

In the simulations, we have fixed some parameters and varied others as
follows. The length and width of the unit cell are fixed as $l_{x}=12\xi$ and
$l_{y}=8\xi$. The thickness of the sample varied from $l_{z}=1\xi$ to
$l_{z}=3.6\xi$ in increments of $0.2\xi$. We use
$\kappa=1/\sqrt{2},1,\sqrt{3}$ for each set of values of $l_{z}$. The grid
space used is $\Delta x=\Delta y=0.2\xi$ and $\Delta z=0.1\xi$. The size of
the simulation box was chosen sufficiently large in order to satisfy boundary
conditions (9); we use $l_{y}=16\xi$ and $l_{z}=12\xi$. The dimensions of the
defects are $a_{x}=a_{y}=0.2\xi$. The range $10\leq\gamma\leq 20$ is suitable
for most metals like Nb [17, 18, 23]; we used $\gamma=10$.

Let us explain how we calculate the IV (current-voltage) and IR (current-
resistance) characteristics, which are the measurable quantities in the
resistive state. In the Weyl gauge, the electrical field is given by
$\textbf{E}=-\partial\textbf{A}/\partial t$, so, assuming that the voltage is
measured between electrodes at $z=0$ that cover the width of the film, the
voltage across a unit cell is:

$\displaystyle U(t)$ $\displaystyle=$
$\displaystyle-\frac{1}{n_{y}-1}\sum_{j=2}^{n_{y}}\int_{-l_{x}/2}^{l_{x}/2}\,E_{x}(x,y_{j},0)dx$
(10) $\displaystyle=$
$\displaystyle\frac{1}{n_{y}-1}\sum_{j=2}^{n_{y}}\int_{-l_{x}/2}^{l_{x}/2}\frac{\partial
A_{x}(x,y_{j},0)}{\partial t}dx,$

where $n_{y}=l_{y}/\Delta y$, and $y_{j}=(j-n_{y}/2-1)\Delta y$ for are the
$y$ coordinates of the mesh points. The voltage is then calculated as a time
average of $U(t)$. We have:

$V=\frac{1}{\cal T}\,\int_{0}^{\cal T}\,U(t)\,dt,$ (11)

where ${\cal T}$ is the time corresponding to an appropriate number of
oscillations of $U(t)$.

The applied current density was adiabatically increased in steps of $\Delta
J_{a}=0.01J_{GL}$ from the Meissner state until the superconductivity was
fully destroyed. In the resistive state, we moved from value of $J_{a}$ to
$J_{a}+\Delta J_{a}$ only after the voltage $U(t)$ became periodic, which is
the same periodicity with which the v-av pairs are formed and annihilated.
When multiple nucleations of v-av (vortex-antivortex) pairs are present, the
voltage looses its periodicity, and therefore we change the value of $J_{a}$
only after 220 oscillations of $U(t)$ in order to obtain a more accurate value
of the time average voltage.

The results of all simulations are compiled in the following Subsections.

###  III.2 Field Profile, Closed Vortex, and Current Distribution

Figure 2: The magnetic field profile in the vertical plane $x=0$ (parallel to
the $yz$ plane): for better visualization purposes, the arrows are not in real
size; the rectangle inside is a cross section of the superconductor; this
picture is for $\kappa=1$, $l_{y}=8\xi$, $l_{z}=3\xi$; the value of the
current density is $J_{a}=0.26J_{GL}$ just before the critical current density
$J_{c1}=0.27J_{GL}$. The vortex (antivortex) nucleates on the right-hand side
(left hand-side) of the figure.

In Fig. 2 we illustrate the vector field profile at $x=0$ plane (in the middle
of the unit cell, where the defects are located). The local field has symmetry
as if the current density was uniform. Due to the geometric symmetry and the
demagnetization effect, the local magnetic field is larger near the surface of
the superconductor, and decreases deep inside. As can also be observed, the
field is larger on the lateral sides of the sample. In this figure, we show
the vector field for a value of the current density just before the first
critical current $J_{c1}=0.27J_{GL}$. Therefore, once the resistive state sets
in, it is on the lateral sides that the vortex and the antivortex sprout, move
to the center, and finally annihilate each other at the center of the sample.
Then, a periodic collapsing of v-av pairs is established.

Figure 3: Color maps of the superconducting Cooper-pair density, $|\psi|$, for
$\kappa=1$, $l_{z}=3\xi$, and $J_{a}=0.27J_{GL}$ throughout the plane $x=0$:
(a) a v-av pair of curved vortices; on the left hand-side (right hand side) is
the antivortex (vortex); (b) a combination of a vortex and an antivortex
producing a closed vortex; (c) a closed vortex diminishing its radius; (d) a
closed vortex shrinking down at the center. The dark strips on both sides are
due to the defects. These pictures correspond to the same region highlighted
in Fig. 2.

Next, we discuss the morphology of a closed vortex in the resistive state. The
first works about closed vortices were conducted on long current-carrying
superconducting cylinders, so that the vortex follows the geometry of the
sample since from the surface until it collapses at the center [4, 7]. In the
present scenario, we deal with a film of rectangular cross section. Thus,
before the closed vortex is formed, two curved vortices (a vortex and an
antivortex) nucleate in opposite sides of the sample (see panel (a) of Fig. 3)
and move toward the center. Then, as they encounter each other, their ends
join together forming a closed vortex (panel (b)). Once this ringlike vortex
is formed, its radius starts decreasing (panel (c)) until it collapses at the
center. After the transition from the Meissner state to the resistive one, the
process is repeated periodically until superconductivity is suppressed
throughout the sample.

Figure 4: The same as in Fig.3 for $l_{z}=1.6\xi$ and $J_{a}=0.30J_{GL}$.

It is conceivable that the closed vortex can exist for any thickness of the
film, but for small values of $l_{z}$ is very elongated when the ends of the
v-av pair meet. In our simulations, we do not have sufficient resolution to
detect a closed vortex for any $l_{z}$. Indeed, in Fig. 4 we show four panels
of the color maps of the superconducting Cooper-pair density for
$l_{z}=1.6\xi$. As we can see, the shape of the closed vortex is much more
elongated than for the previous case $l_{z}=3\xi$ of Fig. 3. For $\kappa=1$,
and thickness below $l_{z}=1.6\xi$, we do not observe any closed vortex; the
v-av pair remains straight lines, since the nucleation of the vortex and the
antivortex on the surfaces, until the pair is annihilated at the center of the
sample. The formation of a closed vortex depends on the sample thickness
because the current and the magnetic field concentrate near the surface of the
superconductor, on a scale of the order of the London penetration depth.
Therefore, for thicker samples, the ends of a vortex line are subjected to
stronger Lorentz force and to a stronger horizontal magnetic field that
facilitate the formation of a closed vortex.

Figure 5: The panels show four cuts of the current distribution of the closed
vortex throughout the vertical plane $y=0$. The radius of the closed vortex
diminishes from (a) to (d).

Let us now discuss the current distribution of a closed vortex. When both
curved vortex and antivortex touch their ends on the upper and lower surfaces,
$z=+l_{z}/2$ and $z=-l_{z}/2$, respectively, they combine in order to make a
single closed vortex. This new vortex looks like a toroid with the
superconducting currents flowing around its core. Fig. 5 exhibits four cuts of
the toroid in the $xz$ plane for $l_{z}=3\xi$. As it can be seen, the currents
in the internal parts of the toroid flow in the same direction for both the
upper and lower segment of the closed vortex. Therefore, all v-av pairs
opposedly positioned in the toroid attract one another, causing the closed
vortex to collapse at the center of the sample.

###  III.3 Straight to Curved Vortex Crossover, and (IV,IR) Characteristics

Figure 6: (a) Voltage across the $z$ direction as a function of the thickness
of the sample for three values of $\kappa$: the value of $J_{a}$ for each case
corresponds to the first critical current density when the resistive state
sets in. The highlighted dots are the critical $l_{z,c}$ values for which the
v-av pairs combine to make a closed vortex. (b) The derivative of the voltage:
the dots separate the two regimes of straight to curved vortices; the
inflection points coincide with $l_{z,c}$.

As we mentioned previously, as the thickness of the sample is increased, there
is a crossover between straight to curved v-av pair. In what follows, we show
a consistence between the criterion based on the aspect ratio of the vortex
(antivortex) and an important physical quantity, namely, the voltage across
the $z$ direction on the lateral side of the film; by aspect ratio, we mean
the distance between the center of the vortex and the antivortex along the $y$
direction when their tips first touch each other. For this purpose, we
calculate the time average of the following voltage:

$\displaystyle U(t)$ $\displaystyle=$
$\displaystyle-\frac{2}{n_{x}}\sum_{i=n_{x}/2}^{n_{x}}\int_{0}^{l_{z}/2}\,E_{z}(x_{i},l_{y}/2,z)dz$
(12) $\displaystyle=$
$\displaystyle\frac{2}{n_{x}}\sum_{i=n_{x}/2}^{n_{x}}\int_{0}^{l_{z}/2}\frac{\partial
A_{z}(x_{i},l_{y}/2,z)}{\partial t}dz,$

where $n_{x}=l_{x}/\Delta x$, and $x_{i}=(i-n_{x}/2-1)\Delta x$ for all
$\\{i=1,2,\ldots,n_{x}+1\\}$ are the $x$ coordinates of the mesh points. Here,
we have not considered the branch $-l_{z}/2\leq z\leq 0$. By symmetry, had we
included this contribution, the voltage would vanish.

When a closed vortex appears, we will have a larger contribution for the
current flowing in the vertical direction, and consequently an increase in the
voltage. For a fixed value of $\kappa$ and current density $J_{a}$, we
determine the voltage for several values of $l_{z}$. We have done this for
three distinct values of the Ginzburg-Landau parameter. The respective value
of $J_{a}$ is chosen so as to correspond to the critical current density
$J_{c1}$ for the lowest thickness, $l_{z}=1\xi$. The results are summarized in
Fig. 6. In panel (a) we present the voltage as a function of the thickness of
the film. We find that, at a certain point, which we denote by $l_{z,c}$,
there is a change of the behavior of the $V(l_{z})$ curves. These points are
highlighted in panel (a). They signal a crossover from linear to curved v-av
pairs.

In order to make sure that this special point is correlated to the straight-
to-curved vortex crossover, we calculate the derivative $dV(l_{z})/dl_{z}$ for
the three values of $\kappa$ (see panel (b)). As can be seen, the derivatives
have an inflection point which are highlighted in panel (b). These points
correspond to $l_{z,c}$. We find the following critical values,
$l_{z,c}=1.4\xi,1.6\xi,1.8\xi$ for $\kappa=1/\sqrt{2},1,\sqrt{3}$,
respectively. We must emphasize that, first we determine the value of
$l_{z,c}$ by inspecting the aspect ratio of the curvature of the vortex.
Second, we check if the result is in agreement with inflection point of
$dV(l_{z})/dl_{z}$. For all the three cases mentioned above they coincide.

Figure 7: IV (blue line) and IR (red line) characteristic curves,
respectively, for $\kappa=1$ and $l_{z}=3\xi$. The Meissner state (full
superconductivity) survives up to $J_{a}=J_{c1}=0.27J_{GL}$. Above this
current density, the resistive state sets in. The resistive state splits into
two phases. In one of them the vortex and the antivortex nucleate only at the
defects on the border of the superconductor. In the second phase, another set
of v-av pairs nucleates at the frontiers between unit cells. The second jump
in the IV characteristic is the signature of this crossover. The insets
illustrate this scenario through the modulus of the order parameter in the
$xy$ plane ($z=0$ plane).

Now we discuss the transport properties of the superconductor. The IV and IR
characteristics are presented in Fig. 7 for $\kappa=1$ and $l_{z}=3\xi$. As we
increase the applied current density, the system becomes unstable to the
penetration of v-av pairs. When $J_{a}$ achieves the value $J_{c1}=0.27J_{GL}$
the superconductor goes to the resistive state, where a periodic formation
v-av pairs occurs. Notice that $J_{c1}$ is smaller than the depairing current
density $J_{GL}$. This is a consequence of the defects deliberately introduced
at the border of the superconducting film. If we further increase $J_{a}$, a
second jump appears in the IV curve at $J_{a}=0.575J_{GL}$. This is an
indication that another two adjacent v-av pairs around the central one are
nucleating (see insets). This is in correspondence with the experimental
observations of multiple penetrations of kinematic vortices in Sn film by
Sivakov et al. [24]. Finally, when the current density reaches the value
$J_{a}=J_{c2}=1.155J_{GL}$ the superconductor goes straight to the normal
state. We believe that for larger unit cells, additional jumps in the IV curve
would occur.

###  III.4 Single Defect (Half-Closed Vortex)

Figure 8: Color maps of the superconducting Cooper-pair density, $|\psi|$, for
$\kappa=1$, $l_{z}=4\xi$, and $J_{a}=0.26J_{GL}$ throughout the plane $x=0$:
(a) an av nucleates in the left edge of the sample and moves towards the
opposite side; (b) the ends of the av touch the $y=l_{y}/2$ plane and form a
half-closed vortex; (c) a half-closed vortex diminishing its radius; (d) the
half-closed vortex shrinking down.

We also have considered a single defect in the middle of a unit cell (in the
middle of front edge in Fig. 1). In this configuration, only an antivortex
nucleates on the $y=-l_{y}/2$ surface. Once the antivortex nucleates at
$y=-l_{y}/2$ it moves directly towards the opposite side. As can be seen from
Fig. 8, as the antivortex approaches the other side of the sample, it becomes
significantly curved (see panel (a)). When it reaches the surface $y=l_{y}/2$,
surprisingly, it does not escape the sample. Instead, its ends touch the
surface giving rise to a half-closed vortex (see panels (b) and (c)). Then, it
diminishes its ratio until it collapses (see panel (d)).

Due to its intrinsic nature, it is very difficult to observe experimentally
the closed vortex. In addition, both the closed and half-closed vortex are
very unstable. Therefore, we require an indirect method that signals either a
v-av or a single half-closed vortex curves. Having this in mind, we propose a
setup to detect the curvature of the v-av pair when it gives rise to a closed
vortex with a non vanishing aspect ratio. Instead of doing this for a closed
vortex, which is formed inside the sample, we think it should be much easier
for a half-closed vortex, since its collapse occurs on the surface. For this
purpose, we calculate the time average of the of the magnetic flux on the
lateral side of the film. In order to calculate the magnetic flux, we focus in
a small region where the antivortex tips touch the plane $y=l_{y}/2$, although
we could extend it throughout the whole lateral side of the unit cell. We
evaluate the following equation:

$\Phi(t)=\int_{0}^{l_{z}/2}\int_{-\xi}^{\xi}\,h_{y}(x,\pm l_{y}/2,z)\,dxdz.$
(13)

Here, the minus sign stands for the left edge where the antivortex nucleates,
and the plus sign is for the opposite one where the antivortex ends touch the
surface. We consider only half of the lateral edge, otherwise the total flux
would vanish.

Figure 9: The main panel presents the time average of the magnetic flux,
$\bar{\Phi}$, across the surfaces defined in Eq. 13; the black points
correspond to the beginning of the resistive state. The inset shows the
magnetic flux, $\Phi$, as a function of time. The parameters used were
$\kappa=1$, $l_{z}=4\xi$, and $J_{a}=0.26J_{GL}$.

Fig. 9 presents the results for the time averaged magnetic flux by using the
same parameters as those used in Fig. 8. As can be clearly seen, until the
transition to the resistive state, the flux is approximately the same through
both surfaces $y=\pm l_{y}/2$. Nevertheless, once the resistive state sets in,
they become different as much as $\Delta\bar{\Phi}\approx
0.09\xi^{2}H_{c2}\approx 0.02\Phi_{0}$. This is signaling that the antivortex
is piercing the $y=+l_{y}/2$ surface. Therefore, for this to happen, the
antivortex necessarily has to bend.

Since the creation and annihilation of the half-closed vortex is a dynamical
process, the flux evolves periodically. The AC magnetic flux can be seen in
the inset of Fig. 9. The period of the AC signal depends on the applied
current density. For $J_{a}=0.26J_{GL}$ we find that the period is
$\tau\approx 0.03\times 10^{4}t_{GL}$. For low-$T_{c}$ materials like Nb films
[23], $t_{GL}\approx 6.72$ ps. This produces $\tau\approx 2$ ns, which is in
the GHz frequency range.

As we can see, the measurement of the difference between the time averaged
magnetic flux threading at each plane, as displayed in Fig. 9, can be an
indirect method for the experimental detection of a closed vortex. Such
measurement is experimentally feasible by using the recently developed
nanoSQUIDs [25, 26, 27], which are capable of detecting the variation of the
flux produced by the closed vortex in the time and length scale we used in our
computations.

Figure 10: Difference between the average magnetic flux on both sides of the
sample through two vertical circuits positioned on the planes
$y=\pm(l_{y}/2+\xi)$. The domain of the circuits is given by $\\{-\xi\leq
x\leq\xi,\,-l_{z}/2\leq z\leq L_{z}/2\\}$. The value of $L_{z}$ was chosen
such that the area of the circuit above the $z=l_{z}/2$ surface is the same
for all thicknesses $l_{z}$. The points just before the onset of the resistive
state are highlighted in black.

We must emphasize that for thin superconductors, the measurement of the flux
in the region prescribed in the above setup can be experimentally challenging.
For this reason, we also present another indirect method for the detection of
a closed vortex. Fig 10 shows the difference between the averaged magnetic
flux calculated at the lateral sides of our superconductor as a function of
the applied current density, for different superconductor thicknesses. In
contrast with the previous case, here the flux is calculated from the bottom
of the superconductor up to a height well above the sample surface. The flux
is evaluated across a vertical rectangular surface, located a coherence length
away from the lateral surface. This makes the proposed experiment much more
feasible.

The black dots in Fig. 10 represent the onset of the resistive state for each
thickness. We are interested in current densities slightly above these values.
In this region, the antivortex moves through the whole sample, being expelled
at the other side (animations of this regime, as well as the one described
below, can be found in the Supplement Material 111See Supplemental Material at
[URL] for animations of the dynamics of the half-closed vortex.). Due to its
curvature, the antivortex produces a larger magnetic flux in the plane which
it is moving into, increasing the flux difference between each plane. Since
the curvature increases with the film thickness, this difference also
increases with the sample thickness, as shown in Fig. 10. Nevertheless, for
large values of the current density, a vortex also penetrates the
superconductor at the opposite side, with the pair being annihilated inside
the sample, reducing the impact of the curvature in the flux. The penetration
of this vortex becomes easier as the superconductor thickness increases, which
explains why the flux difference for $l_{z}=4\xi$ is smaller than for
$l_{z}=2\xi$ or $3\xi$ at high current densities, for example. In summary, by
comparing the magnetic flux difference, at the onset of the resistive state,
for films with different thicknesses, we can clearly demonstrate the existence
of the half-closed vortex. Given the inherent complexity for the direct
observation of a closed vortex, our indirect method brings a new possibility
for the first detection of such objects.

##  IV Concluding remarks

To summarize, we have shown that the combination of the flux lines of a vortex
and an antivortex during their annihilation gives origin to a closed vortex
loop. As we show here, the formation of the closed vortex depends on how
easily the flux lines can be bent due to the action of the applied current,
with this bending increasing with the superconducting film thickness and
decreasing with the Ginzburg-Landau parameter $\kappa$. Since the motion and
annihilation of vortices are highly dissipative processes, understanding their
behavior is of fundamental importance to the design of electronic devices.

Our findings suggest a new method to experimentally observe a closed vortex.
As discussed here, closed vortices can be indirectly detected by the measuring
the flux produced by their stray fields. We emphasize that the recently
developed nanoSQUIDs are capable of performing such measurements in the time
and length scales that our system requires.

###### Acknowledgements.

LRC, LVT, and ES thank the Brazilian Agency FAPESP for financial support,
grant numbers 20/03947-2, 19/24618-0, 20/10058-0, respectively. ES thanks
Professor Alexey Samokhvalov for very useful discussions. ES is also grateful
for the warming hospitality of the Department of Physics, University of
Antwerp, where this work was finished. WAO thanks the National Council for
Scientific and Technological Development (CNPq, Grant 309928/2018-4). ES
thanks Professor Felipe Fernandes Fanchini for kindly donating two GPU cards
(FAPESP, grant number 21/04655-8).

## References

  * Abrikosov [1957] A. A. Abrikosov, On the Magnetic properties of superconductors of the second group, Sov. Phys. JETP 5, 1174 (1957), [Zh. Eksp. Teor. Fiz. 32,1442(1957)].
  * Kozlov and Samokhvalov [1991] V. Kozlov and A. Samokhvalov, Closed Abrikosov vortices in type-II superconductors, Pis’ma Zh. Eksp. Teor. Fiz. 53 (1991).
  * Kozlov and Samokhvalov [1993a] V. Kozlov and A. Samokhvalov, Closed Abrikosov vortices in a superconducting cylinder, Physica C: Superconductivity and its Applications 213, 103 (1993a).
  * Kozlov and Samokhvalov [1993b] V. Kozlov and A. Samokhvalov, Stabilization of toroidal Abrikosov vortex in a nonuniform superconductor, Journal of Superconductivity 6 (1993b).
  * Genenko [1993] Y. Genenko, Relaxation of magnetic vortex rings in a superconducting cylinder: some universal features, Physica C: Superconductivity 215, 343 (1993).
  * Genenko [1994a] Y. A. Genenko, Vortex helicoid in a superconducting cylinder in a longitudinal magnetic field, JETP Letters C/C of Pis’ma Zh. Eksp. Teor. Fiz. 59, 841 (1994a).
  * Genenko [1994b] Y. A. Genenko, Magnetic self-field entry into a current-carrying type-ii superconductor, Phys. Rev. B 49, 6950 (1994b).
  * Genenko [1995] Y. A. Genenko, Magnetic self-field entry into a current-carrying type-II superconductor. II. Helical vortices in a longitudinal magnetic field, Phys. Rev. B 51, 3686 (1995).
  * Genenko  _et al._ [1998] Y. A. Genenko, A. V. Snezhko, P. Troche, J. Hoffmann, and H. C. Freyhardt, Magnetic self-field entry into a current-carrying type-II superconductor. III. general criterion of penetration for an external field of arbitrary direction, Phys. Rev. B 57, 1164 (1998).
  * Samokhvalov [1996] A. Samokhvalov, Vortex loops entry into type-II superconductors, Physica C: Superconductivity 259, 337 (1996).
  * Samokhvalov [1997] A. Samokhvalov, Abrikosov vortex loop near the surface of superconductor, Physica C: Superconductivity and its Applications 282-287, 2163 (1997).
  * Samokhvalov [1998] A. V. Samokhvalov, Expanding vortex rings in a current-carrying superconducting cylinder, Physica C: Superconductivity 308, 74 (1998).
  * Niedzielski and Berakdar [2022] B. Niedzielski and J. Berakdar, Vortex ring and helical current formation in superconductors driven by a thz-field-induced toroidal vector potential, physica status solidi (b) , 2100622 (2022).
  * Berdiyorov  _et al._ [2018] G. Berdiyorov, M. Milošević, F. Kusmartsev, F. Peeters, and S. Savel’ev, Josephson vortex loops in nanostructured josephson junctions, Scientific reports 8, 1 (2018).
  * Fyhn and Linder [2019] E. H. Fyhn and J. Linder, Controllable vortex loops in superconducting proximity systems, Physical Review B 100, 214503 (2019).
  * Gurevich [2017] A. Gurevich, Theory of rf superconductivity for resonant cavities, Superconductor Science and Technology 30, 034004 (2017).
  * Kramer and Watts-Tobin [1978] L. Kramer and R. Watts-Tobin, Theory of dissipative current-carrying states in superconducting filaments, Physical Review Letters 40, 1041 (1978).
  * Watts-Tobin  _et al._ [1981] R. Watts-Tobin, Y. Krähenbühl, and L. Kramer, Nonequilibrium theory of dirty, current-carrying superconductors: Phase-slip oscillators in narrow filaments near ${T}_{c}$, Journal of Low Temperature Physics 42, 459 (1981).
  * Du and Gray [1996] Q. Du and P. Gray, High-$\kappa$ limits of the time-dependent Ginzburg-Landau model, SIAM Journal on Applied Mathematics 56, 1060 (1996).
  * Cadorim  _et al._ [2020] L. R. Cadorim, A. de Oliveira Junior, and E. Sardella, Ultra-fast kinematic vortices in mesoscopic superconductors: the effect of the self-field, Scientific Reports 10, 1 (2020).
  * Barba-Ortega  _et al._ [2015] J. Barba-Ortega, E. Sardella, and J. A. Aguiar, Superconducting properties of a parallelepiped mesoscopic superconductor: A comparative study between the 2D and 3D Ginzburg–Landau models, Physics Letters A 379, 732 (2015).
  * Gropp  _et al._ [1996] W. D. Gropp, H. G. Kaper, G. K. Leaf, D. M. Levine, M. Palumbo, and V. M. Vinokur, Numerical simulation of vortex dynamics in type-II superconductors, Journal of Computational Physics 123, 254 (1996).
  * Berdiyorov  _et al._ [2009] G. R. Berdiyorov, M. V. Milošević, and F. M. Peeters, Kinematic vortex-antivortex lines in strongly driven superconducting stripes, Physical Review B 79 (2009).
  * Sivakov  _et al._ [2003] A. G. Sivakov, A. M. Glukhov, A. N. Omelyanchouk, Y. Koval, P. Müller, and A. V. Ustinov, Josephson behavior of phase-slip lines in wide superconducting strips, Phys. Rev. Lett. 91, 267001 (2003).
  * Embon  _et al._ [2017] L. Embon, Y. Anahory, Ž. L. Jelić, E. O. Lachman, Y. Myasoedov, M. E. Huber, G. P. Mikitik, A. V. Silhanek, M. V. Milošević, A. Gurevich, _et al._ , Imaging of super-fast dynamics and flow instabilities of superconducting vortices, Nature communications 8, 1 (2017).
  * Anahory  _et al._ [2020] Y. Anahory, H. Naren, E. Lachman, S. B. Sinai, A. Uri, L. Embon, E. Yaakobi, Y. Myasoedov, M. Huber, R. Klajn, _et al._ , Squid-on-tip with single-electron spin sensitivity for high-field and ultra-low temperature nanomagnetic imaging, Nanoscale 12, 3174 (2020).
  * Vasyukov  _et al._ [2013] D. Vasyukov, Y. Anahory, L. Embon, D. Halbertal, J. Cuppens, L. Neeman, A. Finkler, Y. Segev, Y. Myasoedov, M. L. Rappaport, _et al._ , A scanning superconducting quantum interference device with single electron spin sensitivity, Nature nanotechnology 8, 639 (2013).
  * Note [1] See Supplemental Material at [URL] for animations of the dynamics of the half-closed vortex.


# Strategies to reduce the thermoelastic loss of multimaterial coated finite
substrates

R. Zhou ID Department of Materials Science and Engineering, University of
California, Berkeley, California 94720, USA  M. Molina-Ruiz ID Department of
Physics, University of California, Berkeley, California 94720, USA  F. Hellman
ID Department of Materials Science and Engineering, University of California,
Berkeley, California 94720, USA  Department of Physics, University of
California, Berkeley, California 94720, USA

(August 27, 2024)

###### Abstract

Thermoelastic loss is one of the main energy dissipation mechanisms in
resonant systems. A careful analysis of the thermoelastic loss is critical to
design low-noise resonators for high-precision applications, such as
gravitational-wave detectors. This paper presents an analytical solution to
the thermoelastic loss in multimaterial coated finite substrates with
realistic assumptions on the model structure and the elastic fields. The
mechanism responsible for thermoelastic loss is taken as a function of
material properties, operating temperature and frequency, and other design
parameters. We calculate the thermoelastic loss for specific applications over
a wide range of frequencies ($1$ Hz to 10 GHz) and temperatures ($1$ K to
$300$ K), and for a variety of substrate and coating materials. The result is
relevant for gravitational-wave detectors and for experiments sensitive to
mechanical dissipation.

###### pacs:

##  I Introduction

Thermoelastic dissipation arises from the thermal response generated by the
elastic field in a deformable medium. In a vibrating structure, the periodic
displacement field causes a temperature gradient in the elastic body, where
the compressed region is hotter and the stretched region becomes cooler.
Consequently, heat transfer takes place in order to reach thermal equilibrium,
which results in elastic components out of phase with the elastic fields. The
energy lost due to this irreversible heat flow is referred to as thermoelastic
loss [1, 2]. For a coated substrate, this loss is also caused by thermal
expansion mismatch between the materials that coating and substrate are made
of. This type of coating-substrate structure is commonly used in mirrors used
for the detection of gravitational-waves (GW), such as in the Laser
Interferometer Gravitational-Wave Observatory (LIGO) [3], Virgo interferometer
[4] and the Kamioka Gravitational-Wave Detector (KAGRA) [5], for which it is
necessary to estimate the loss generated by this effect.

Thermoelastic loss of gravitational-wave mirrors has been extensively studied
in addition to other sources of thermal noise, such as the mechanical loss
intrinsic to the material and thermo-refractive noise caused by temperature
fluctuations [6]. The quality factor, Q-factor, is a measure of the energy
loss rate and is the inverse measure of energy loss caused by all the
dissipation mechanisms present in the structure. To further improve the
Q-factor of resonant systems, or the ratio of stored energy to the energy
dissipated per cycle, it is essential to fully understand the thermoelastic
loss and predict it for various experimental conditions.

In 1937, Zener first published an analysis of thermoelastic damping in thin
resonator beams undergoing flexural vibrations [1]. Lifshitz and Roukes later
refined his model and solved the fundamental equations more rigorously,
leading to different approximate solutions of temperature profile and hence,
different expressions of the thermoelastic loss (hereinafter referred to as
$\phi$) [7]. After solving the classical Fourier heat conduction equation,
they derived $\phi$ from the attenuation of the vibration.

Fejer et al. developed an independent approach for the computation of $\phi$
in mirrors for gravitational-wave detectors; they determined the energy lost
due to the elastic field induced by thermal fluctuations in the coating on an
infinitely thick substrate [8]. Based on the work of Fejer et al., Somiya and
Yamamoto considered a substrate of finite thickness and proposed a different
form of solution to the thermal equations using the elastic response of a
cylinder with finite thickness [9]. The energy dissipation resulting from the
heat flow was directly extracted from the temperature gradient.

Gravitational-wave detectors, such as LIGO, are kilometer-sized
interferometers that bounce a beam of light between highly reflective optical
mirrors. These mirrors consist of a silica substrate and a multilayer coating
that alternates layers of silica and titania-doped tantala [10]. It is thus of
high relevance to calculate the thermoelastic loss in mirrors with these types
of structures. Fejer et al. and Somiya and Yamamoto, both adopted the same
approach to calculate the thermoelastic loss in mirrors with multilayered
coatings: the effective medium approach (EMA), in which the multilayers are
abstracted as a homogeneous medium of weighted-average physical properties [8,
9].

In this paper, we present a theoretical treatment of the thermoelastic
behavior of a coated substrate, based on the model described by Fejer et al.
[8]. We propose that a more realistic determination of $phi$ can be achieved
using a model that i) considers a substrate of finite thickness, and ii)
preserves the system configuration and the identity of the materials that
compose it. The derived expressions are then applied to various structures and
as will be shown in this paper, the thermoelastic loss can be minimized
through careful choices of variables, such as operating frequency,
temperature, and thickness of individual coating layers. It should be noted
that our model includes only the thermoelastic loss resulting from the
interaction between two or more materials of different physical properties,
when they are in contact with each other; the internal loss within one
material due to statistical fluctuations in temperature has been studied for
finite-sized mirrors and can be combined with the coupled loss discussed in
this paper to find the total thermoelastic loss [11, 12, 13].

The paper is outlined as follows: Section II presents an overview of the key
equations used to evaluate the thermoelastic loss in a substrate coated by a
uniform film. Thermoelastic response and the consequent heat propagation
normal to the surfaces are studied. We also consider a model consisting of a
multilayer coating and refine the heat equations to accommodate the
alternating layers. Analytical solutions for both types of structure are
derived. Finally, Section III presents the results calculated from the
formulae of thermoelastic loss for various cases and discusses how material
properties and external factors would impact the thermoelastic loss.
Particular attention is paid to the multilayer structure without adopting the
EMA, where the structure of the coating (layers and interfaces) is preserved.

##  II Mathematical framework for the thermoelastic loss

To derive an expression for the thermoelastic loss in a coated substrate, we
consider a film of thickness $l$ on a substrate with a thickness of $h-l$, as
illustrated in Figure 1. We define the surface normal to be in the
$z$-direction. The surface is located at $z=0$, and the film extends from
$z=0$ to $z=l$, and the substrate from $z=l$ to $z=h$, where $h\gg l$. Our
model assumes that both the film and the substrate are uniform so that there
is no temperature variation in the $x$-$y$ plane, and the transverse
dimensions are much larger than the thickness so that only thermal diffusion
along $z$-direction needs to be taken into account.

Figure 1: Illustration of a uniform substrate coated with a layer of
homogeneous material.

The loss in energy due to thermoelasticity of materials is calculated for two
states of stress, in-plane stress and normal stress. Detailed derivations are
presented in the following sections. Thermoelastic loss is a measure of energy
dissipated due to the coupling between the deformation and thermal field.
Hence first, the relationship between the intrinsic elastic field and the
temperature field produced by it needs to determined. By defining the boundary
conditions and solving the equilibrium equations, the thermal strain and the
associated power dissipation can be found. Subsequently, the elastic energy
stored in the coating can be obtained from the applied stress through the use
of theory of elasticity. Finally, based on the definition of thermoelastic
loss $\phi$, an expression of $\phi$ for the coating and the substrate can be
found.

###  II.1 Formulation of the thermal field

The dynamic deformation of a body causes temperature variations. We use the
linear heat equation along $z$-direction to find the temperature distribution
coupled to the input elastic field [14],

$\frac{\partial\theta_{j}}{\partial
t}-\kappa_{j}\frac{\partial^{2}\theta_{j}}{\partial
z^{2}}=-\frac{E_{j}\alpha_{j}T}{(1-2\nu_{j})C_{j}}\frac{\partial}{\partial
t}\sum\epsilon_{0,j}$ (1)

where $\theta_{j}(t,z)$ is the time/position-varying temperature, $\kappa_{j}$
is the thermal diffusitivity, $E_{j}$ is the Young’s modulus, $\alpha_{j}$ is
the coefficient of linear thermal expansion, $T$ is the background
temperature, $\nu_{j}$ is the Poisson’s ratio, $C_{j}$ is the specific heat
capacity per unit volume and $j=f$, $s$ indicates quantities evaluated in the
film and in the substrate, respectively. GW mirrors are exposed to high-
intensity laser beams, which causes that the mirror thermal field depends on
its elastic field and optical absorption. In the model presented in this
paper, the thermal field only depends on the elastic field. The work on
photothermal transfer function by Ballmer [15] serves as a basis for
calculating the heat flow and thermal expansion caused by optical absorption.

Taking $\theta_{j}(t,z)$ to be in the form of $\theta_{j}(z)exp(i\omega t)$
and the strain as $\epsilon_{0}exp(i\omega t)$, where $\omega=2\pi f$ is the
angular frequency, Equation 1 becomes,

$i\omega\theta_{j}-\kappa_{j}\frac{\partial^{2}\theta_{j}}{\partial
z^{2}}=-i\omega\beta_{j}$ (2)

where $\beta_{j}=(E_{j}\alpha_{j}T)/[(1-2\nu_{j})C_{j}]\sum\epsilon_{0,j}$.

The boundary conditions for heat fluxes are defined as follows,

$\displaystyle\frac{\partial\theta_{f}}{\partial
z}\Bigr{|}_{\begin{subarray}{c}z=0\end{subarray}}$ $\displaystyle=0\text{ ,}$
(3) $\displaystyle k_{f}\frac{\partial\theta_{f}}{\partial
z}\Bigr{|}_{\begin{subarray}{c}z=l\end{subarray}}$
$\displaystyle=k_{s}\frac{\partial\theta_{s}}{\partial
z}\Bigr{|}_{\begin{subarray}{c}z=l\end{subarray}}\text{ ,}$
$\displaystyle\frac{\partial\theta_{s}}{\partial
z}\Bigr{|}_{\begin{subarray}{c}z=h\end{subarray}}$ $\displaystyle=0$

where $k$ represents the thermal conductivity.

As previously proposed by Fejer et al. [8], the final solution to the heat
equation consists of two parts, a particular (p) solution of thermal field,
$\theta_{p,j}(z)$, that is coupled to the strain and satisfies Equation 2, and
a specific (s) solution, $\theta_{s,j}(z)$, that meets the boundary
conditions. It can be written as,

$\theta_{j}=\theta_{p,j}(z)+\theta_{s,j}(z)$ (4)

We solve for the thermal fields consistent with the boundary conditions
(Equation 3); the calculation is shown in Appendix A. The thermal fields in
the film and the substrate along the $z$ direction are

$\displaystyle\theta_{f}$ $\displaystyle=-\beta_{f}+\frac{\Delta\beta\times
cosh(\gamma_{f}z)}{\cosh(\gamma_{f}l)+R\sinh(\gamma_{f}l)\coth(q)}\text{ ,}$
(5) $\displaystyle\theta_{s}$ $\displaystyle=-\beta_{s}-\frac{\Delta\beta
R\times cosh[\gamma_{s}(h-z)]}{\coth(\gamma_{f}l)\sinh(q)+R\cosh(q)}$

where $q=\gamma_{s}(h-l)$, $\Delta\beta=\beta_{f}-\beta_{s}$,
$\gamma_{j}=(1+i)\sqrt{\pi fC_{j}/k_{j}}$ and
$R=(k_{f}\gamma_{f})/(k_{s}\gamma_{s})=\sqrt{(k_{f}C_{f})/(k_{s}C_{s})}$.

To determine the induced stress field, $\epsilon_{1}$, and strain field,
$\sigma_{1}$, by oscillatory thermal fluctuation, we assume that any
expansion/contraction along $z$-direction is not restricted, and therefore,
$\epsilon_{1,zz}=\alpha_{j}\theta_{j}$, $\sigma_{1,zz}=0$. When considering
the in-plane thermal expansion, one assumption made in the calculation is that
the substrate expands freely, effectively uninfluenced by the film; however,
the in-plane expansion of the film is constrained. Moreover, the film is
assumed to have a uniform strain of $\alpha_{s}\theta_{s}(z=l)$, which is the
expansion of the substrate at the film-substrate interface. This approximation
is reasonable if the substrate is much thicker than the film, i.e., $(h-l)\gg
l$. Since the free expansion in plane of the film is forbidden, stress is
developed in the interior of the film and the stress level depends on the
mismatch between the thermal expansion of the film and the substrate. Taking
these into consideration, we find the elastic fields due to thermal mismatch:

$\displaystyle\epsilon_{1,ii,j}$
$\displaystyle=A_{1,ii,j}\alpha_{j}\theta_{j}\text{ ,}$ (6)
$\displaystyle\sigma_{1,ii,j}$ $\displaystyle=B_{1,ii,j}\alpha_{j}\theta_{j}$

The matrices $A_{1,ii,j}$ and $B_{1,ii,j}$ are defined as,

$\displaystyle A_{1,ii,f}$
$\displaystyle=\begin{pmatrix}\alpha_{s}\theta_{s,l}&0&0\\\
0&\alpha_{s}\theta_{s,l}&0\\\ 0&0&\alpha_{f}\theta_{f}\end{pmatrix}\text{ ,}$
(7) $\displaystyle B_{1,ii,f}$ $\displaystyle=\begin{pmatrix}-b&0&0\\\
0&-b&0\\\ 0&0&0\end{pmatrix}\text{ ,}$ $\displaystyle A_{1,ii,s}$
$\displaystyle=\begin{pmatrix}\alpha_{s}\theta_{s}&0&0\\\
0&\alpha_{s}\theta_{s}&0\\\ 0&0&\alpha_{s}\theta_{s}\end{pmatrix}\text{ ,}$
$\displaystyle B_{1,ii,s}$ $\displaystyle=\begin{pmatrix}0&0&0\\\ 0&0&0\\\
0&0&0\end{pmatrix}$

where $b=E_{f}(\alpha_{f}\theta_{f}-\alpha_{s}\theta_{s,l})$, and
$\theta_{s,l}$ represents $\theta_{s}(z=l)$.

###  II.2 Applied elastic field and elastic energy

In order to calculate the oscillatory thermal field caused by an oscillatory
applied stress or a vibration at some frequency $\omega$, we need to first
determine the stress and strain states. Here two potential cases are
considered: 1) in-plane stress and 2) normal stress to the surface. The latter
one is relevant to mechanical loss measurements and LIGO operation. Solutions
to most other stress states can be attained by a sum of the solutions based on
these two stress fields. Using the elastic boundary conditions, we find the
driving elastic fields in the film, denoted by the subscript $0$.

####  II.2.1 Stress parallel to the coated surface

The boundary conditions for stress parallel to the coated surface are:
$\sigma_{0,xx}=\sigma_{0,yy}=\sigma_{0,\parallel}$ (with the assumption of
symmetry), and $\sigma_{0,zz}=0$. Using Hooke’s law, we find

$\displaystyle\epsilon_{0,xx,j}=$
$\displaystyle\epsilon_{0,yy,j}=\frac{1-v_{j}}{E_{j}}\sigma_{0,\parallel}\text{
,}$ (8) $\displaystyle\epsilon_{0,zz,j}$
$\displaystyle=\frac{-2v_{j}}{E_{j}}\sigma_{0,\parallel}$

Summarizing the results in matrix form, the stress and strain fields can be
expressed in the forms of

$\displaystyle\epsilon_{0,ii,\parallel,j}$
$\displaystyle=a_{0,ii,j}\sigma_{0,\parallel}\text{ ,}$ (9)
$\displaystyle\sigma_{0,ii,\parallel,j}$
$\displaystyle=b_{0,ii,j}\sigma_{0,\parallel}$

where

$\displaystyle a_{0,j}$
$\displaystyle=\begin{pmatrix}\frac{1-v_{j}}{E_{j}}&0&0\\\
0&\frac{1-v_{j}}{E_{j}}&0\\\ 0&0&\frac{-2v_{j}}{E_{j}}\end{pmatrix}\text{ ,}$
(10) $\displaystyle b_{0,j}$ $\displaystyle=\begin{pmatrix}1&0&0\\\ 0&1&0\\\
0&0&0\end{pmatrix}$

Therefore the elastic energy stored per unit area is given by

$\displaystyle E_{stored,\parallel,j}$
$\displaystyle=\frac{1}{2}\sigma\epsilon\times l$ (11)
$\displaystyle=\frac{1}{2}\sigma_{0,\parallel,j}^{2}l\sum b_{0,j}a_{0,j}$
$\displaystyle=\frac{1}{2}\sigma_{0,\parallel,j}^{2}l\times\frac{2(1-v_{j})}{E_{j}}$

####  II.2.2 Stress perpendicular to the coated surface

For stress perpendicular to the coated surface, the boundary conditions are
defined as: $\sigma_{0,zz}=\sigma_{0,\perp}$, $\sigma_{0,xx},\sigma_{0,yy}\neq
0$; $\epsilon_{0,zz}\neq 0$, $\epsilon_{0,xx}=\epsilon_{0,yy}=0$. With the
Hooke’s law expression, we obtain

$\displaystyle\sigma_{0,xx,j}$
$\displaystyle=\sigma_{0,yy,j}=\frac{\nu_{j}}{1-\nu_{j}}\sigma_{0,\perp}\text{
,}$ (12) $\displaystyle\epsilon_{0,zz,j}$
$\displaystyle=\frac{\sigma_{0,\perp}}{E_{j}}\frac{(1-2\nu_{j})(1+\nu_{j})}{1-\nu_{j}}$

The collective results are,

$\displaystyle\epsilon_{0,ii,\perp,j}$
$\displaystyle=c_{0,ii,j}\sigma_{0,\perp}\text{ ,}$ (13)
$\displaystyle\sigma_{0,ii,\perp,j}$
$\displaystyle=d_{0,ii,j}\sigma_{0,\perp}$

where

$\displaystyle c_{0,j}$ $\displaystyle=\begin{pmatrix}0&0&0\\\ 0&0&0\\\
0&0&\frac{(1-2\nu_{j})(1+\nu_{j})}{1-\nu_{j}}\frac{1}{E_{j}}\end{pmatrix}\text{
,}$ (14) $\displaystyle d_{0,j}$
$\displaystyle=\begin{pmatrix}\frac{\nu_{j}}{1-\nu_{j}}&0&0\\\
0&\frac{\nu_{j}}{1-\nu_{j}}&0\\\ 0&0&1\end{pmatrix}$

Hence, the amount of elastic energy stored per unit area is given by

$\displaystyle E_{stored,\perp,j}$
$\displaystyle=\frac{1}{2}\sigma\epsilon\times l$ (15)
$\displaystyle=\frac{1}{2}\sigma_{0,\perp,j}^{2}l\sum d_{0,j}c_{0,j}$
$\displaystyle=\frac{1}{2}\sigma_{0,\perp,j}^{2}l\times\frac{(1-2\nu_{j})(1+\nu_{j})}{E_{j}(1-\nu_{j})}$

###  II.3 Energy dissipation and thermoelastic loss

We first consider the energy dissipated within the film. The rate of energy
dissipation per unit volume in a deformed body, in this case, the film, is
defined as

$p_{diss,f}=\frac{1}{2}\frac{\sigma\partial\epsilon}{\partial t}$ (16)

which has to be a real function. The overall film stress and strain in the
case of plane stress are described as

$\displaystyle\epsilon_{f,ii,\parallel}$
$\displaystyle=a_{0,ii}\epsilon_{0}+A_{1,ii,f}\alpha_{f}\theta_{f}\text{ ,}$
(17) $\displaystyle\sigma_{f,ii,\parallel}$
$\displaystyle=b_{0,ii}\sigma_{0}+B_{1,ii,f}\alpha_{f}\theta_{f}$

As for the case of stress perpendicular to the coated surface, we write

$\displaystyle\epsilon_{f,ii,\perp}$
$\displaystyle=c_{0,ii}\epsilon_{0}+A_{1,ii,f}\alpha_{f}\theta_{f}\text{ ,}$
(18) $\displaystyle\sigma_{f,ii,\perp}$
$\displaystyle=d_{0,ii}\sigma_{0}+B_{1,ii,f}\alpha_{f}\theta_{f}$

The oscillatory stress and strain fields have a time dependence

$\displaystyle\epsilon=\epsilon_{f,ii}\times e^{i\omega t}\text{ ,}$ (19)
$\displaystyle\sigma=\sigma_{f,ii}\times e^{i\omega t}$

Hence, $p_{diss,f}$ can be written as

$\displaystyle p_{diss,f}$
$\displaystyle=\sum\frac{1}{2}\frac{\partial\sigma_{f,ii}\epsilon_{f,ii}}{\partial
t}$ (20)
$\displaystyle=\omega\sum\operatorname{Im}[\sigma_{f,ii}\epsilon_{f,ii}]$
$\displaystyle=\omega\sum\operatorname{Im}[\sigma_{0,ii}\epsilon_{0,ii}+\sigma_{1,ii}\epsilon_{0,ii}$
$\displaystyle\quad+\sigma_{0,ii}\epsilon_{1,ii}+\sigma_{1,ii}\epsilon_{1,ii}]$
$\displaystyle=\omega\sum\operatorname{Im}[\sigma_{1,ii}\epsilon_{0,ii}+\sigma_{0,ii}\epsilon_{1,ii}]$

since $\sigma_{0,ii}\epsilon_{0,ii}$ is real, and
$\sigma_{1,ii}\epsilon_{1,ii}$ can be neglected because the induced elastic
fields are significantly smaller than the input fields.

Therefore, the energy dissipated per unit area is

$E_{diss,f}=\tau\int_{0}^{l}p_{diss,f}dz$ (21)

where $\tau=2\pi/\omega$. Substituting Equation 20 and the expressions of
thermal and applied elastic fields, previously derived in Section II.1, into
Equation 21, we obtain the total energy dissipated in the film and the
substrate for parallel and perpendicular stress fields, respectively:

$\displaystyle E_{diss,total,\parallel}$
$\displaystyle=2\pi\sigma_{0}\times\operatorname{Im}\left[(2\nu_{f}-2)\alpha_{f}\theta_{1,f}\gamma_{f}^{-1}\sinh(\gamma_{f}l)\right.$
(22)
$\displaystyle\quad+(4-2\nu_{f})\alpha_{s}\theta_{1,s}\cosh(\gamma_{s}(h-l))\times
l$
$\displaystyle\quad-2\alpha_{s}\theta_{1,s}\gamma_{s}^{-1}\sinh(\gamma_{s}(h-l))\Big{.}\Big{]}$
$\displaystyle E_{diss,total,\perp}$
$\displaystyle=2\pi\sigma_{0}\times\operatorname{Im}\left[\alpha_{f}\theta_{1,f}\gamma_{f}^{-1}\sinh(\gamma_{f}l)\right.$
(23)
$\displaystyle\quad+\frac{2\nu_{f}}{1-\nu_{f}}\alpha_{s}\theta_{1,s}\cosh[\gamma_{s}(h-l)]\times
l$
$\displaystyle\quad-\frac{1+\nu_{s}}{1-\nu_{s}}\alpha_{s}\theta_{1,s}\gamma_{s}^{-1}\sinh[\gamma_{s}(h-l)]\Big{.}\Big{]}$

The detailed derivation is included in Appendix B.

The loss factor is generally defined as the ratio of energy dissipated per
radian to the potential energy in a cycle. In this case, the potential energy
is the elastic energy stored in the strained material. The consideration of an
infinite substrate prevents the analysis of the elastic energy stored in the
substrate, and only the strain energy of the coating can be considered.
Therefore, to allow a direct addition/subtraction of all loss components
(mechanical and thermoelastic), the denominator has to be consistently defined
for these loss factors as the elastic energy of the entire system (coating and
substrate). In this paper, we consider the stored potential energy in both the
coating and the substrate.

According to Equations 11 and 15, the stored energy in the entire system can
be written in the form of

$\displaystyle E_{stored,total,\parallel}$
$\displaystyle=\frac{1}{2}\sigma_{0}^{2}l\times\frac{2(1-\nu_{f})}{E_{f}}$
(24)
$\displaystyle\quad+\frac{1}{2}\sigma_{0}^{2}(h-l)\times\frac{2(1-\nu_{s})}{E_{s}}\text{
,}$ $\displaystyle E_{stored,total,\perp}$
$\displaystyle=\frac{1}{2}\sigma_{0}^{2}l\times\frac{(1-2\nu_{f})(1+\nu_{f})}{E_{f}(1-\nu_{f})}$
$\displaystyle\quad+\frac{1}{2}\sigma_{0}^{2}(h-l)\times\frac{(1-2\nu_{s})(1+\nu_{s})}{E_{s}(1-\nu_{s})}$

Therefore, the thermoelastic loss $\phi$ can thus be calculated from

$\displaystyle\phi_{\parallel}$
$\displaystyle=\frac{\left|E_{diss,total,\parallel}\right|}{2\pi\times
E_{stored,total,\parallel}}\text{ ,}$ (25) $\displaystyle\phi_{\perp}$
$\displaystyle=\frac{\left|E_{diss,total,\perp}\right|}{2\pi\times
E_{stored,total,\perp}}$

Typically the substrate is significantly thicker than the film, i.e.,
$(h-l)\gg l$, implying that the total stored energy is mostly contained in the
substrate (the second terms of Equations 24). As a result, by using a thicker
substrate, which leads to a larger elastic energy stored in the substrate
(while having little impact on the temperature gradient and inducing only a
small increase in the energy lost per cycle due to thermoelasticity), we see a
decrease in $\phi$. It should be noted that this is different from the
dilution factor considered in mechanical loss extraction from GeNS measurement
[16], since the thermoelastic loss is coupled, due to the interactions between
coating and substrate, and cannot be separated. Somiya and Yamamoto showed
that $\phi$ decreases with the substrate thickness and hits a plateau at a
certain thickness [9]. As per the definition of thermoelastic loss (Eq. 25),
adding more material to the substrate increases the system stored energy,
while the energy dissipated remains almost unchanged. Thus, using a thicker
substrate would lead to a lower thermoelastic loss. The thermoelastic loss
plateau discussed above is a direct consequence of only taking into account
the film elastic energy. The substrate effect will be discussed in more detail
in Section III.

###  II.4 Multilayer coating on a finite substrate

In this section, we consider a coating made of alternating layers of two
materials attached on a homogeneous substrate, which resembles the dielectric
mirror coating used in gravitational-wave detectors, as depicted in Figure 2.
We make the assumptions that i) the outer surface is not subjected to any heat
flux, namely, that heat transfer only takes place within the films and
substrate, and ii) each component is only affected by its nearest neighbors.
These assumptions are likely to be valid in most cases as the temperature
fluctuation caused by the driving stress is reasonably small so that there
would not exist a large temperature difference across the coated body.
Therefore, the thermoelastic dissipation can be divided into two parts: energy
loss due to thermal conduction between 1) the coating layers, and 2) between
the first deposited layer and the substrate.

Focusing on the interaction in the multilayer coating, we assume that each
layer is uniform and the system is in steady state. To find the energy
dissipation between coating layers made of two different materials, we begin
by considering the thermal field in a single layer and assuming that its
derivative is zero at the center, a direct consequence of continuity. The
thermal fields are then solved using the boundary conditions and a model
consisting of two half-layers [depicted in Figure 2(a)], where the temperature
profile and heat flux are deemed symmetric.

Figure 2: Illustration of a multilayer coating on a uniform substrate: (a)
depicts a schematics of the model used to calculate the loss between layers,
and (b) to calculate the loss between the first deposited layer and the
substrate, where $N$ represents the number of layers in the coating, $t_{1}$
the thickness of a single coating layer made of Material 1, $t_{2}$ the
thickness of a single coating layer made of Material 2, and $t_{s}$ the
substrate thickness. The red lines show an example of the temperature profile
across coating and substrate.

The thermal field solutions are similar to those derived in Section II.1, with
the boundary conditions now being

$\displaystyle\frac{\partial\theta_{f_{1}}}{\partial
z}\Bigr{|}_{\begin{subarray}{c}z=0\end{subarray}}$ $\displaystyle=0\text{ ,}$
(26) $\displaystyle k_{f_{1}}\frac{\partial\theta_{f_{1}}}{\partial
z}\Bigr{|}_{\begin{subarray}{c}z=\frac{t_{1}}{2}\end{subarray}}$
$\displaystyle=k_{f_{2}}\frac{\partial\theta_{f_{2}}}{\partial
z}\Bigr{|}_{\begin{subarray}{c}z=\frac{t_{1}}{2}\end{subarray}}\text{ ,}$
$\displaystyle\frac{\partial\theta_{f_{2}}}{\partial
z}\Bigr{|}_{\begin{subarray}{c}z=\frac{t_{1}}{2}+\frac{t_{2}}{2}\end{subarray}}$
$\displaystyle=0$

Solving for the thermal fields in the two materials, we get

$\displaystyle\theta_{f_{1}}$
$\displaystyle=\theta_{p,f_{1}}+\theta_{1,f_{1}}cosh(\gamma_{f_{1}}z)\text{
,}$ (27) $\displaystyle\theta_{f_{2}}$
$\displaystyle=\theta_{p,f_{2}}+\theta_{1,f_{2}}cosh\left[\gamma_{f_{2}}(\frac{t_{1}}{2}+\frac{t_{2}}{2}-z)\right]$

And $\theta_{1,f_{1}}$ and $\theta_{1,f_{2}}$ are defined as

$\displaystyle\theta_{1,f_{1}}$
$\displaystyle=\frac{\Delta\beta^{{}^{\prime}}}{\cosh(\frac{\gamma_{f_{1}}t_{1}}{2})+R^{{}^{\prime}}\sinh(\frac{\gamma_{f_{1}}t_{1}}{2})\coth(\frac{\gamma_{f_{2}}t_{2}}{2})}\text{
,}$ (28) $\displaystyle\theta_{1,f_{2}}$
$\displaystyle=-\frac{\Delta\beta^{{}^{\prime}}R^{{}^{\prime}}}{\coth(\frac{\gamma_{f_{1}}t_{1}}{2})\sinh(\frac{\gamma_{f_{2}}t_{2}}{2})+R^{{}^{\prime}}\cosh(\frac{\gamma_{f_{2}}t_{2}}{2})}$

where $\Delta\beta^{{}^{\prime}}=\beta_{f_{1}}-\beta_{f_{2}}$ and
$R^{{}^{\prime}}=(k_{f_{1}}\gamma_{f_{1}})/(k_{f_{2}}\gamma_{f_{2}})=\sqrt{(k_{f_{1}}C_{f_{1}})/(k_{f_{2}}C_{f_{2}})}$.
The energy dissipated for parallel and perpendicular fields in the adjacent
films (illustrated in Figure 2(a), indicated by subscript $a$, respectively,
is thus obtained:

$\displaystyle E_{diss,a,\parallel}$
$\displaystyle=2\pi\sigma_{0}\times\operatorname{Im}\left[(2\nu_{f1}-2)\alpha_{f1}\theta_{1,f1}\gamma_{f1}^{-1}\sinh\left(\frac{\gamma_{f1}t_{1}}{2}\right)\right.$
(29)
$\displaystyle\quad+(4-2\nu_{f1})\alpha_{s}\theta_{1,s}\cosh(\gamma_{s}t_{s})\times\frac{t_{1}}{2}$
$\displaystyle\quad-(2-2\nu_{f2})\alpha_{f2}\theta_{1,f2}\gamma_{f2}^{-1}\sinh\left(\frac{\gamma_{f2}t_{2}}{2}\right)$
$\displaystyle\quad+(4-2\nu_{f2})\alpha_{s}\theta_{1,s}\cosh(\gamma_{s}t_{s})\times\frac{t_{2}}{2}\bigg{.}\bigg{]}$
$\displaystyle E_{diss,a,\perp}$
$\displaystyle=2\pi\sigma_{0}\times\operatorname{Im}\left[\alpha_{f1}\theta_{1,f1}\gamma_{f1}^{-1}\sinh\left(\frac{\gamma_{f1}t_{1}}{2}\right)\right.$
(30)
$\displaystyle\quad+\frac{2\nu_{f1}}{1-\nu_{f1}}\alpha_{s}\theta_{1,s}\cosh(\gamma_{s}t_{s})\times\frac{t_{1}}{2}$
$\displaystyle\quad+\alpha_{f2}\theta_{1,f2}\gamma_{f2}^{-1}\sinh\left(\frac{\gamma_{f2}t_{2}}{2}\right)$
$\displaystyle\quad+\frac{2\nu_{f2}}{1-\nu_{f2}}\alpha_{s}\theta_{1,s}\cosh(\gamma_{s}t_{s})\times\frac{t_{2}}{2}\bigg{.}\bigg{]}$

From elasticity theory, the elastic energy stored in the case of parallel and
perpendicular fields, respectively, in these half layers can be expressed as

$\displaystyle E_{stored,a,\parallel}$
$\displaystyle=\frac{1}{2}\sigma_{0}^{2}\frac{t_{1}}{2}\times\frac{2(1-\nu_{f1})}{E_{f1}}$
(31)
$\displaystyle\quad+\frac{1}{2}\sigma_{0}^{2}\frac{t_{2}}{2}\times\frac{2(1-\nu_{f2})}{E_{f2}}\text{
,}$ $\displaystyle E_{stored,a,\perp}$
$\displaystyle=\frac{1}{2}\sigma_{0}^{2}\frac{t_{1}}{2}\times\frac{(1-2\nu_{f_{1}})(1+\nu_{f_{1}})}{E_{f_{1}}(1-\nu_{f_{1}})}$
$\displaystyle\quad+\frac{1}{2}\sigma_{0}^{2}\frac{t_{2}}{2}\times\frac{(1-2\nu_{f_{2}})(1+\nu_{f_{2}})}{E_{f_{2}}(1-\nu_{f_{2}})}$

Similarly, for the thermoelastic loss due to heat conduction between the last
layer of Material $2$ and the substrate (illustrated in Figure 2(b), indicated
by subscript $b$, we can find the energy dissipated in the form of

$\displaystyle E_{diss,b,\parallel}$
$\displaystyle=2\pi\sigma_{0}\times\operatorname{Im}\left[(2\nu_{f2}-2)\right.$
(32)
$\displaystyle\quad\times\alpha_{f2}\theta^{{}^{\prime}}_{1,f2}\gamma_{f2}^{-1}\sinh\left(\frac{\gamma_{f2}t_{2}}{2}\right)$
$\displaystyle\quad+(4-2\nu_{f2})\alpha_{s}\theta_{1,s}\cosh(\gamma_{s}t_{s})\times\frac{t_{2}}{2}$
$\displaystyle\quad-2\alpha_{s}\theta_{1,s}\gamma_{s}^{-1}\sinh(\gamma_{s}t_{s})\bigg{.}\bigg{]}$
$\displaystyle E_{diss,b,\perp}$
$\displaystyle=2\pi\sigma_{0}\times\operatorname{Im}\left[\alpha_{f2}\theta^{{}^{\prime}}_{1,f2}\gamma_{f2}^{-1}\sinh\left(\frac{\gamma_{f2}t_{2}}{2}\right)\right.$
(33)
$\displaystyle\quad+\frac{2\nu_{f2}}{1-\nu_{f2}}\alpha_{s}\theta_{1,s}\cosh(\gamma_{s}t_{s})\times\frac{t_{2}}{2}$
$\displaystyle\quad-\frac{1+\nu_{s}}{1-\nu_{s}}\alpha_{s}\theta_{1,s}\gamma_{s}^{-1}\sinh(\gamma_{s}t_{s})\bigg{.}\bigg{]}$

where $\theta^{{}^{\prime}}_{1,f_{2}}$ and $\theta_{1,s}$ are defined as

$\displaystyle\theta^{{}^{\prime}}_{1,f_{2}}$
$\displaystyle=\frac{\Delta\beta^{{}^{\prime\prime}}}{\cosh(\frac{\gamma_{f_{2}}t_{2}}{2})+R^{{}^{\prime\prime}}\sinh(\frac{\gamma_{f_{2}}t_{2}}{2})\coth(\gamma_{s}t_{s})}\text{
,}$ (34) $\displaystyle\theta_{1,s}$
$\displaystyle=-\frac{\Delta\beta^{{}^{\prime\prime}}R^{{}^{\prime\prime}}}{\coth(\frac{\gamma_{f_{2}}t_{2}}{2})\sinh(\gamma_{s}t_{s})+R^{{}^{\prime\prime}}\cosh({\gamma_{s}t_{s})}}$

and where $\Delta\beta^{{}^{\prime\prime}}=\beta_{f_{2}}-\beta_{s}$ and
$R^{{}^{\prime\prime}}=(k_{f_{2}}\gamma_{f_{2}})/(k_{s}\gamma_{s})=\sqrt{(k_{f_{2}}C_{f_{2}})/(k_{s}C_{s})}$.

And the energy stored for parallel and perpendicular fields in the form of

$\displaystyle E_{stored,b,\parallel}$
$\displaystyle=\frac{1}{2}\sigma_{0}^{2}\frac{t_{2}}{2}\times\frac{2(1-\nu_{f_{2}})}{E_{f_{2}}}$
(35)
$\displaystyle\quad+\frac{1}{2}\sigma_{0}^{2}t_{s}\times\frac{2(1-\nu_{s})}{E_{s}}$
$\displaystyle E_{stored,b,\perp}$
$\displaystyle=\frac{1}{2}\sigma_{0}^{2}\frac{t_{2}}{2}\times\frac{(1-2\nu_{f_{2}})(1+\nu_{f_{2}})}{E_{f_{2}}(1-\nu_{f_{2}})}$
$\displaystyle\quad+\frac{1}{2}\sigma_{0}^{2}t_{s}\times\frac{(1-2\nu_{s})(1+\nu_{s})}{E_{s}(1-\nu_{s})}$

The $E_{diss}$ and $E_{stored}$ in the first deposited layer of Material $2$
are double counted in the two models. Taking it into account, the final
expression of the total thermoelastic loss of a multilayer coating on a
uniform substrate is given by

$\phi=\frac{E_{diss,a}\times
N-E_{diss,a_{2}}+E_{diss,b}}{2\pi\times(E_{stored,a}\times
N-E_{stored,a_{2}}+E_{stored,b})}$ (36)

where $N$ is the total number of layers, $E_{diss,a_{2}}$ and
$E_{stored,a_{2}}$ stand for the energy dissipated and stored in half a layer
of Material 2, correspondingly.

##  III Results and Discussion

In this section, we present estimates of thermoelastic loss for various cases,
which are calculated using the equations presented in previous sections. The
input elastic field, unless specified, uses the normal stress field described
in Section II.2, as it is commonly seen in many applications and used in the
measurement of coating properties. The physical properties of materials used
in the calculation are listed in Table 1.

Table 1: Elastic and thermal properties at room temperature of materials used
for the calculation of $\phi$, Young’s modulus $E$, Poisson’s ratio $\nu$,
thermal expansion coefficient $\alpha$, specific heat per unit volume $C_{V}$,
and thermal conductivity $k$. The temperature dependencies of the thermal
properties ($\alpha$, $C_{V}$ and $k$) are implemented as reported in their
corresponding references.

Materials |  $E$ |  $\nu$ |  $\alpha$ |  $C_{V}$ |  $k$
---|---|---|---|---|---
|  (GPa) |  |  (10-6 K-1) |  (106 J m-3 K-1) |  (W m-1 K-1)
PVD $a$-Si |  115 a |  0.18 a |  2.4 b [17] |  2.4 [18] |  91.0 [19]
c-Si |  169 [20] |  0.28 [20] |  2.6 [21] |  1.6 [22] |  92.0 [23]
PVD $a$-SiO2 |  66 [24] |  0.17 [24] |  0.6 c [24] |  2.2 c [24] |  1.2 c [24]
$a$-Ti:Ta2O5 |  143 [25] |  0.23 [26] |  3.9 [25] |  2.1 [8] |  33 [8]
$c$-GaAs |  85.5 [27] |  0.31 [27] |  5.75 [27] |  1.74 [27] |  46 [27]
$c$-Al0.92Ga0.08As |  100 [28] |  0.32 [28] |  5.2 [28] |  1.7 [28] |  70 [28]

  * a

PVD $a$-Si elastic properties are obtained using $E=2G(1+\nu)$, where
$E=(1-\nu)140$ GPa [29], and the shear modulus $G=49$ GPa [30], both
experimentally measured.

  * b

PVD $a$-Si thermal expansion coefficient temperature dependence is assumed to
be the same than for c-Si reported in Ref. [21].

  * c

PVD $a$-SiO2 thermal properties temperature dependence is assumed to be the
same as for fused $a$-SiO2 reported in Refs. [31, 32].

###  III.1 Thickness and frequency dependence

The thermoelastic loss of single-layer coated substrates was first
investigated. $\phi$ of amorphous silicon ($a$-Si) films of various
thicknesses, deposited on fused silica and crystalline Si (c-Si) substrates at
room temperature (RT) is plotted in Figure 3 as a function of frequency.
Numerical values of $\phi$ can be obtained from Equation 25, for a frequency
range from $1$ to $10^{7}$ Hz. The calculated $\phi$ shows a non-monotonic
behavior with one or two inflection points as a function of frequency, with
the peak shifting towards lower frequency as the film thickness increases.

Figure 3: Thermoelastic loss as a function of frequency calculated at room
temperature for an $a$-Si coating of different thicknesses on different
substrates that are 100 $\upmu$m thick: (a) fused silica (a-SiO2) and (b)
crystalline silicon (c-Si).

For the effect of thickness, our results show that a thicker film would yield
higher thermoelastic loss, which is a direct outcome of a larger temperature
difference, and thus a higher heat flux between the film and the substrate.
The substrate thickness also has an impact on $\phi$: thicker substrates lead
to an approximately proportional reduction in the thermoelastic loss, which is
demonstrated in Figure 4. As discussed at the end of Section II.3, while the
elastic energy stored is linearly proportional to the substrate thickness,
there is only a marginally small increase in the overall energy dissipation,
which results in near proportional damping. Eventually, $\phi$ would approach
zero at infinite substrate thickness, as the energy loss caused by thermal
expansion mismatch would be vanishingly small in comparison to the elastic
energy stored in the substrate.

Figure 4: Thermoelastic loss as a function of frequency calculated at room
temperature for a 100 nm-thick $a$-Si coating onto silica substrates 100
$\upmu$m-thick (black line) and 100 cm-thick (red line).

Considering that both film and substrate thicknesses have an impact on the
thermoelastic loss, we note that certain film to substrate thickness ratios
may cause the thermoelastic loss to be one of the main energy dissipation
mechanisms in the system.

###  III.2 Effect of input stress field

As previously mentioned in Section II.2, another common type of stress is in-
plane stress. Figure 5 compares the thermoelastic loss for the in-plane and
perpendicular stress components. It can be seen that the peak frequency and
the peak loss are only slightly affected; the difference is within $3\%$ and
$11\%$, respectively. Other types of bulk stress and shear stress can also be
incorporated into the calculation, which makes the model presented in this
work versatile.

Figure 5: Thermoelastic loss as a function of frequency calculated at room
temperature for $a$-Si s with different thicknesses on $100$ $\upmu$m-thick
c-Si substrates with normal stress (solid lines) and in-plane stress (dashed
lines).

###  III.3 Multilayer coating

We now consider a multilayer coating consisting of amorphous silicon and
amorphous silica on top of a fused silica substrate. The substrate thickness
is set to be $100$ $\mu$m and the layer thicknesses are selected to follow the
quarter-wavelength rule, $\lambda/4n$, where $\lambda$ is the wavelength of
the operating laser and $n$ is the refractive index [33]. We first
investigated the effect of the number of layers. As can be seen in Figure 6,
while the peak frequency remains almost unaffected, the thermoelastic loss
increases in proportion to the number of layers with constant layer thickness,
or in other words, the amount of material.

Figure 6: Thermoelastic loss as a function of frequency calculated at room
temperature for a multilayer coating that consists of alternating layers of
PVD $a$-Si (142.9 nm) and $a$-SiO2 (347.2 nm) on a fused silica substrate. We
used $n=3.5$ for $a$-Si [34] and $n=1.44$ for $a$-SiO2 [35]

The thermoelastic loss models from Fejer et al. [8] and Somiya and Yamamoto
[9] consider the EMA to approximate the two materials used in the multilayer
coating as one homogeneous film with their averaged physical properties. In
contrast, our method treats each layer as an individual component and computes
the loss generated from heat transport across all interfaces. In Figure 7 we
compare the results obtained by the model presented in this work and the EMA.
The total film thickness was kept constant while the layer thicknesses were
varied. It can be seen that, the three peaks of the loss curves are
approximately at the same level, which is as expected since $\phi$ is
proportional to the total volume of the system. However, our model predicts
that the number of interfaces plays a role and affects the peak position in
frequency; as the layer thickness decreases, or equivalently, the number of
interfaces increases, the peak position is shifted towards higher frequencies,
which effectively lowers $\phi$ in the frequency range of interest for GW
detectors (100 Hz to 10 kHz), such as LIGO, Virgo and KAGRA. These results
suggest that by reducing the layers thickness without changing the overall
thickness of the coating, increasing the number of interfaces, the
thermoelastic loss at a certain frequency can be reduced.

Figure 7: Thermoelastic loss as a function of frequency calculated at room
temperature for a $100$ $\upmu$m fused silica substrate coated with
alternating $a$-Si and silica layers with a total thickness of $4901$ nm. We
show the results for one layer using the EMA (blue line), two layers where the
$a$-Si layer is 1429 nm-thick and the $a$-SiO2 layer is 3472 nm-thick (black
line), 20 layers where the $a$-Si layers are 142.9 nm-thick and the $a$-SiO2
layers are 347.2 nm-thick (red line), and 200 layers where the $a$-Si layers
are 14.29 nm-thick and the $a$-SiO2 layers are 34.72 nm-thick (orange line).

###  III.4 Comparison of Mirror Coatings for Gravitational-wave Detectors

Based on the model presented in Section II.4, the thermoelastic loss of
multilayer coating is calculated for various combinations of materials and
layer thicknesses at three excitation frequencies, 100 Hz, 1 kHz and 10 kHz.
The predicted values are listed in Table 2. Our calculation shows that the
choice of mirror coating and substrate in Advanced LIGO (the current LIGO
generation); using a coating made of amorphous titania-doped tantala and
amorphous silica and a fused silica substrate, has the lowest thermoelastic
loss at all three frequencies compared to other candidate materials.

Table 2: Estimated thermoelastic loss of current Advanced LIGO mirrors and potential candidates. Mirror components | Laser wavelength | $\phi$
---|---|---
(material 1 / material 2 // substrate) |  | 100 Hz | 1 kHz | 10 kHz
$a$-Ti:Ta2O5 / $a$-SiO2 (30 layers) // $a$-SiO2 | 1064 nm | $5.99\times 10^{-15}$ | $5.99\times 10^{-14}$ | $5.99\times 10^{-13}$
$a$-Si / $a$-SiO2 (20 layers) // $c$-Si | 1550 nm | $4.63\times 10^{-11}$ | $3.07\times 10^{-13}$ | $1.40\times 10^{-12}$
2000 nm | $7.71\times 10^{-12}$ | $5.41\times 10^{-13}$ | $2.63\times 10^{-12}$
$c$-GaAs / $c$-AlGaAs (30 layers) // $c$-Si | 1064 nm | $2.79\times 10^{-12}$ | $3.20\times 10^{-13}$ | $1.01\times 10^{-12}$
1550 nm | $6.15\times 10^{-12}$ | $7.02\times 10^{-13}$ | $2.22\times 10^{-12}$
$a$-Ta2O5 / $a$-Si / $a$-SiO2 (30 layers) // $a$-SiO2 [36] | 1064 nm | $6.75\times 10^{-13}$ | $3.01\times 10^{-12}$ | $1.83\times 10^{-11}$

###  III.5 Dependence on Temperature

Thermoelastic loss is closely related to working temperature since the thermal
field induced is affected by T, as seen in Equation 1. In addition, the
thermal properties of materials (coefficient of thermal expansion $\alpha$,
specific heat $C_{V}$, and thermal conductivity $k$) are temperature
dependent. As a result, $\phi$ varies with temperature. The elastic properties
depend much less on temperature than the other properties so here we assume
that they do not change with temperature, and we use measured values where
possible and made estimations when those are not available for the thermal
properties to calculate the temperature dependence of $phi$. We assume that
amorphous silicon has the same thermal expansion coefficient of its
crystalline form [21] and its thermal conductivity and heat capacity are taken
from Refs. [19, 18]. The values for thermal properties of fused silica were
taken from Refs. [37, 31, 32]. Substituting in these temperature-dependent
properties, the relationship between $\phi$ and temperature can then be
analyzed.

Figure 8 illustrates that the thermoelastic loss of $100$ nm-thick $a$-Si
deposited onto a silica substrate of $100$ $\mu$m thick at $10$ K is
significantly lower than that at RT for all frequencies, and the inflection
point shifts towards higher frequencies as temperature decreases.

Figure 8: Thermoelastic loss as a function of frequency calculated at
different temperatures for a layer of $a$-Si $100$ nm-thick onto a silica
substrate of $100$ $\upmu$m thick.

It is known that the thermal expansion coefficient of crystalline silicon,
$\alpha_{c-Si}$, is zero at $17.6$ K and $123.7$ K [21, 38]. Hence, the future
cryogenic gravitational-wave detector, LIGO Voyager, that plans to operate at
$123.7$ K will greatly reduce the thermoelastic loss by making the test mass
and one of the coating materials from silicon, as proposed by Adhikari et al.
[39]. $a$-Si has good mechanical loss and optical reflectivity compared to the
currently used Ta2O5 coating [40, 41, 42]. Future LIGO Voyager plans to use
crystalline silicon for the mirror substrate, and amorphous silicon as the
high-index material and silica as the low-index layer for the mirror coating.
This detector will operate at cryogenic temperatures to further reduce the
loss and improve the sensitivity beyond the detection limits of the current GW
detectors.

We calculated the thermoelastic loss of the proposed multilayer stack of
$a$-Si/$a$-SiO2 films for the LIGO Voyager mirror coating at three different
excitation frequencies. As plotted in Figure 9, the thermoelastic loss
increases with increasing vibration frequency, which agrees with what is
observed in the calculated loss at RT.

Figure 9: Thermoelastic loss as a function of temperature for a stack of 16
alternating layers of $142.9$ nm-thick $a$-Si and $347.2$ nm-thick $a$-SiO2
onto a $10$ cm-thick c-Si substrate, calculated at an excitation frequency of
$100$ Hz (black line), $1$ kHz (red line) and $10$ kHz (blue line). The black
dashed line ($100$ Hz*) has a reversed order of the layers in the coating so
that $a$-Si is in contact with the $c$-Si substrate. The artifact (sharp peak)
seen at at 10 kHz and $\sim 120$ K is due to limited precision of the MATLAB
calculation.

The predicted loss shown in Figure 9 is based on the structure in which
$a$-SiO2 is the last layer and in contact with the substrate. We note that
when $a$-Si is the last layer, the thermoelastic loss is lower at any
temperature; as the discrepancy between the thermal behavior of $a$-Si and
$c$-Si is smaller than that between $a$-SiO2 and $c$-Si (black dashed line in
Figure 9). If the thermal expansion coefficient of $a$-Si is not exactly zero
at the same temperatures as for $c$-Si (17.6 K and 123.7 K), then the $\phi$
two minima seen in Figure 9 would split into two local minima each, one to the
left and one to the right. This could affect the current plans of making 123.7
K the operating temperature for LIGO Voyager.

##  IV Conclusions

We present a mathematical model that estimates the thermoelastic loss
associated with a coating on a substrate of finite thickness. The numerical
results obtained highlight how material properties, measurement temperature
and frequency, and experimental system design affect the thermoelastic loss.
We find that a thicker coating increases the thermoelastic loss, and in
contrast, a thicker substrate reduces it. Additionally, we further extend our
analytical solutions to multilayered structures and demonstrate that
thermoelastic loss correlates with the number of layers (interfaces) for a
given volume, which is a relevant parameter in the design of gravitational-
wave detector mirrors. By increasing the number of layers, while keeping the
total thickness of the system constant, the thermoelastic loss curve shifts
toward higher frequencies. We also calculate and compare the thermoelastic
loss of the current Advanced LIGO mirror coating and potential candidates for
future gravitational-wave detectors. We pay particular attention to the
amorphous silicon/amorphous silica coating designed for LIGO Voyager and
calculate its thermoelastic loss over a range of temperatures and for
different frequencies. We highlight the relevance of the mismatch between the
physical properties of substrate and first layer at their interface.

###### Acknowledgements.

We thank M. M. Fejer for fruitful discussions. We acknowledge the support of
the LSC Center for Coatings Research, jointly funded by the National Science
Foundation and the Gordon and Betty Moore Foundation through Grant No. 6793.

##  Appendix A Thermal field solutions

As discussed in Section II.1, two components are required for the expression
of thermal fields in $z$-direction, the particular solution and the specific
solution. The particular solution is assumed not to depend on in-plane
position if the film/substrate is uniform. Looking at Equation 2, we can
conclude that

$\theta_{p,j}=-\beta_{j}$ (37)

Using the boundary conditions, we can solve for the specific solution

$\displaystyle\theta_{s,f}$
$\displaystyle=\theta_{1,f}\cosh(\gamma_{f}z)\text{ ,}$ (38)
$\displaystyle\theta_{s,s}$ $\displaystyle=\theta_{1,s}\cosh[\gamma_{s}(h-z)]$

and $\gamma_{j}=(1+i)\sqrt{\omega/(2\kappa_{j})}=(1+i)\sqrt{\pi
fC_{j}/k_{j}}$, so that

$\displaystyle\frac{\partial\theta_{p,j}}{\partial
t}-\kappa_{j}\frac{\partial^{2}\theta_{p,j}}{\partial z^{2}}$
$\displaystyle=-i\omega\beta_{j}\text{ ,}$ (39)
$\displaystyle\frac{\partial\theta_{s,j}}{\partial
t}-\kappa_{j}\frac{\partial^{2}\theta_{s,j}}{\partial z^{2}}$
$\displaystyle=0$

The continuity of temperature and thermal flux at $z=l$ requires

$\displaystyle\theta_{p,f}(z=l)+\theta_{1,f}\cosh(\gamma_{f}l)=$ (40)
$\displaystyle\theta_{p,s}(z=l)+\theta_{1,s}\cosh(\gamma_{s}(h-l))\text{ ,}$
$\displaystyle
k_{f}[\theta^{{}^{\prime}}_{p,f}(z=l)+\theta_{1,f}\gamma_{f}\sinh(\gamma_{f}l)]=$
$\displaystyle
k_{s}[\theta^{{}^{\prime}}_{p,s}(z=l)-\theta_{1,s}\gamma_{s}\sinh(\gamma_{s}(h-l))]$

$\theta^{{}^{\prime}}_{p,f}$ and $\theta^{{}^{\prime}}_{p,s}$ are zero for
uniform film and substrate. Substitute in $\theta_{p,f}$ and $\theta_{p,s}$
and solve for $\theta_{1,f}$, we get

$\displaystyle\theta_{1,f}$
$\displaystyle=\frac{\Delta\beta}{\cosh(\gamma_{f}l)+R\sinh(\gamma_{f}l)\coth(q)}\text{
,}$ (41) $\displaystyle\theta_{1,s}$ $\displaystyle=-\frac{\Delta\beta
R}{\coth(\gamma_{f}l)\sinh(q)+R\cosh(q)}$

where $q=\gamma_{s}(h-l)$, $\Delta\beta=\beta_{f}-\beta_{s}$ and
$R=(k_{f}\gamma_{f})/(k_{s}\gamma_{s})=\sqrt{(k_{f}C_{f})/(k_{s}C_{s})}$.

The thermal field solutions can thus be expressed as a sum of the particular
and the specific solutions and be written as

$\displaystyle\theta_{f}$
$\displaystyle=-\beta_{f}+\frac{\Delta\beta\times\cosh(\gamma_{f}z)}{\cosh(\gamma_{f}l)+R\sinh(\gamma_{f}l)\coth(q)}\text{
,}$ (42) $\displaystyle\theta_{s}$ $\displaystyle=-\beta_{s}-\frac{\Delta\beta
R\times\cosh[\gamma_{s}(h-z)]}{\coth(\gamma_{f}l)\sinh(q)+R\cosh(q)}$

##  Appendix B Energy dissipation due to thermoelastic response

Substituting Equation 20 and the expressions obtained previously in Section
II.2 into Equation 21, we get

$\displaystyle E_{diss,f}$
$\displaystyle=\tau\int_{0}^{l}\omega\sum\operatorname{Im}[\sigma_{1,ii}\epsilon_{0,ii}+\sigma_{0,ii}\epsilon_{1,ii}]dz$
(43)
$\displaystyle=\tau\omega\sigma_{0}\int_{0}^{l}\operatorname{Im}[d_{0,f}]dz$

where
$d_{0,f,\parallel}=B_{1}a_{0}+b_{0}A_{1}=-(2-2\nu_{f})\alpha_{f}\theta_{f}+(4-2\nu_{f})\alpha_{s}\theta_{s,l}$
and
$d_{0,f,\perp}=B_{1}c_{0}+d_{0}A_{1}=\alpha_{f}\theta_{f}+2\nu_{f}(\alpha_{s}\theta_{s,l})/(1-\nu_{f})$.

In previous calculations in Section II.1, we have shown that $\theta_{p,f}$ is
a constant and a real number. Hence Equation 43 can be simplified into

$\displaystyle E_{diss,f,\parallel}$
$\displaystyle=2\pi\sigma_{0}\times\operatorname{Im}\left[(2\nu_{f}-2)\frac{\alpha_{f}\theta_{1,f}}{\gamma_{f}}[\sinh(\gamma_{f}l)]\right.$
(44) $\displaystyle\quad+(4-2\nu_{f})\alpha_{s}\theta_{1,s}\cosh(q)\times
l\Big{.}\Big{]}\text{ ,}$ $\displaystyle E_{diss,f,\perp}$
$\displaystyle=2\pi\sigma_{0}\times\operatorname{Im}\left[\frac{\alpha_{f}\theta_{1,f}}{\gamma_{f}}[\sinh(\gamma_{f}l)]\right.$
$\displaystyle\quad+\frac{2\nu_{f}}{1-\nu_{f}}\alpha_{s}\theta_{1,s}\cosh(q)\times
l\Big{.}\Big{]}$

where $q=\gamma_{s}(h-l)$.

A similar derivation for the substrate can be done by following all the
procedures discussed above. The energy dissipated in the substrate can be
determined and expressed in

$\displaystyle E_{diss,s}$
$\displaystyle=\tau\int_{l}^{h}\omega\sum\operatorname{Im}[\sigma_{1,ii,s}\epsilon_{0,ii,s}+\sigma_{0,ii,s}\epsilon_{1,ii,s}]dz$
(45)
$\displaystyle=\tau\omega\sigma_{0}\int_{l}^{h}\operatorname{Im}[d_{0,s}]dz$

where $d_{0,s,\parallel}=2\alpha_{s}\theta_{s}$ and
$d_{0,s,\perp}=(1+\nu_{s})\alpha_{s}\theta_{s}/(1-\nu_{s})$. And

$\displaystyle E_{diss,s,\parallel}$
$\displaystyle=2\pi\sigma_{0}\times\operatorname{Im}\left[-2\frac{\alpha_{s}\theta_{1,s}}{\gamma_{s}}\sinh(q)\right]$
(46) $\displaystyle E_{diss,s,\perp}$
$\displaystyle=2\pi\sigma_{0}\times\operatorname{Im}\left[-\frac{1+\nu_{s}}{1-\nu_{s}}\frac{\alpha_{s}\theta_{1,s}}{\gamma_{s}}\sinh(q)\right]$

In order to calculate the thermoelastic loss in the film and the substrate,
the total energy dissipated and the total elastic energy stored have to be
found. The total energy dissipated is simply the sum of $E_{diss,f}$ and
$E_{diss,s}$ and is given by

$\displaystyle E_{diss,total,\parallel}$
$\displaystyle=E_{diss,f,\parallel}+E_{diss,s,\parallel}$ (47)
$\displaystyle=2\pi\sigma_{0}\times\operatorname{Im}\left[(2\nu_{f}-2)\frac{\alpha_{f}\theta_{1,f}}{\gamma_{f}}\sinh(\gamma_{f}l)\right.$
$\displaystyle\quad+(4-2\nu_{f})\alpha_{s}\theta_{1,s}\cosh(q)\times l$
$\displaystyle\quad-2\frac{\alpha_{s}\theta_{1,s}}{\gamma_{s}}\sinh(q)\bigg{.}\bigg{]}$
$\displaystyle E_{diss,total,\perp}$
$\displaystyle=E_{diss,f,\perp}+E_{diss,s,\perp}$ (48)
$\displaystyle=2\pi\sigma_{0}\times\operatorname{Im}\left[\frac{\alpha_{f}\theta_{1,f}}{\gamma_{f}}(\sinh\gamma_{f}l)\right.$
$\displaystyle\quad+\frac{2\nu_{f}}{1-\nu_{f}}\alpha_{s}\theta_{1,s}\cosh(q)\times
l$
$\displaystyle\quad-\frac{1+\nu_{s}}{1-\nu_{s}}\frac{\alpha_{s}\theta_{1,s}}{\gamma_{s}}\sinh(q)\bigg{.}\bigg{]}$

## References

  * Zener [1937] C. Zener, Physical Review 52, 230 (1937).
  * Zener [1938] C. Zener, Physical Review 53, 90 (1938).
  * Aasi  _et al._ [2015] J. Aasi _et al._ , Classical and Quantum Gravity 32, 074001 (2015).
  * Acernese  _et al._ [2015] F. Acernese _et al._ , Classical and Quantum Gravity 32, 024001 (2015).
  * Aso  _et al._ [2013] Y. Aso, Y. Michimura, K. Somiya, M. Ando, O. Miyakawa, T. Sekiguchi, D. Tatsumi, and H. Yamamoto, Physical Review D 88, 043007 (2013).
  * Evans  _et al._ [2008] M. Evans, S. Ballmer, M. Fejer, P. Fritschel, G. Harry, and G. Ogin, Physical Review D 78, 102003 (2008).
  * Lifshitz and Roukes [2000] R. Lifshitz and M. L. Roukes, Physical Review B 61, 5600 (2000).
  * Fejer  _et al._ [2004] M. M. Fejer, S. Rowan, G. Cagnoli, D. R. M. Crooks, A. Gretarsson, G. M. Harry, J. Hough, S. D. Penn, P. H. Sneddon, and S. P. Vyatchanin, Physical Review D 70, 082003 (2004).
  * Somiya and Yamamoto [2009] K. Somiya and K. Yamamoto, Physical Review D 79, 102004 (2009).
  * Reid and Martin [2016] S. Reid and I. Martin, Coatings 6, 61 (2016).
  * Braginsky  _et al._ [1999] V. B. Braginsky, M. L. Gorodetsky, and S. P. Vyatchanin, Physics Letters, Section A: General, Atomic and Solid State Physics 264, 1 (1999).
  * Liu and Thorne [2000] Y. T. Liu and K. S. Thorne, Physical Review D 62, 122002 (2000).
  * Lovelace  _et al._ [2018] G. Lovelace, N. Demos, and H. Khan, Classical and Quantum Gravity 35, 10.1088/1361-6382/aa9ccc (2018).
  * L.D. Landau [1986] E. L. L.D. Landau,  _Theory of Elasticity_ (Elsevier, 1986).
  * Ballmer [2015] S. W. Ballmer, Physical Review D 91, 023010 (2015).
  * Vajente [2019] G. Vajente,  _Ligo-T1900276_ , Tech. Rep. (2019).
  * De Lima  _et al._ [1999] M. M. De Lima, R. G. Lacerda, J. Vilcarromero, and F. C. Marques, Journal of Applied Physics 86, 4936 (1999).
  * Queen  _et al._ [2013] D. R. Queen, X. Liu, J. Karel, T. H. Metcalf, and F. Hellman, Physical Review Letters 110, 135901 (2013).
  * Zink  _et al._ [2006] B. L. Zink, R. Pietri, and F. Hellman, Physical Review Letters 96, 055902 (2006).
  * Hopcroft  _et al._ [2010] M. A. Hopcroft, W. D. Nix, and T. W. Kenny, Journal of Microelectromechanical Systems 19, 229 (2010).
  * White and Minges [1997] G. K. White and M. L. Minges, International Journal of Thermophysics 18, 1269 (1997).
  * Flubacher  _et al._ [1959] P. Flubacher, A. J. Leadbetter, and J. A. Morrison, Philosophical Magazine 4, 273 (1959).
  * Glassbrenner and Slack [1964] C. J. Glassbrenner and G. A. Slack, Physical Review 134, A1058 (1964).
  * El-Kareh [1995] B. El-Kareh,  _Fundamentals of Semiconductor Processing Technology_ (Springer US, Boston, MA, 1995).
  * Abernathy  _et al._ [2014] M. R. Abernathy, J. Hough, I. W. Martin, S. Rowan, M. Oyen, C. Linn, and J. E. Faller, Applied Optics 53, 3196 (2014).
  * Granata  _et al._ [2016] M. Granata, E. Saracco, N. Morgado, A. Cajgfinger, G. Cagnoli, J. Degallaix, V. Dolique, D. Forest, J. Franc, C. Michel, L. Pinard, and R. Flaminio, Physical Review D 93, 012007 (2016).
  * Sze and Ng [2006] S. Sze and K. K. Ng, in  _Physics of Semiconductor Devices, 3rd Edition_ (John Wiley & Sons, Inc., Hoboken, NJ, USA, 2006) pp. 790–790.
  * Chalermsongsak  _et al._ [2016] T. Chalermsongsak, E. D. Hall, G. D. Cole, D. Follman, F. Seifert, K. Arai, E. K. Gustafson, J. R. Smith, M. Aspelmeyer, and R. X. Adhikari, Metrologia 53, 860 (2016).
  * Witvrouw and Spaepen [1993] A. Witvrouw and F. Spaepen, Journal of Applied Physics 74, 7154 (1993).
  * Molina-Ruiz  _et al._ [2021] M. Molina-Ruiz, Y. J. Rosen, H. C. Jacks, M. R. Abernathy, T. H. Metcalf, X. Liu, J. L. DuBois, and F. Hellman, Physical Review Materials 5, 035601 (2021).
  * Jacobs  _et al._ [1984] S. F. Jacobs, D. Shough, and C. Connors, Applied Optics 23, 4237 (1984).
  * Zeller and Pohl [1971] R. C. Zeller and R. O. Pohl, Physical Review B 4, 2029 (1971).
  * Sheppard [1995] C. J. Sheppard, Pure and Applied Optics: Journal of the European Optical Society Part A 4, 665 (1995).
  * Pierce and Spicer [1972] D. T. Pierce and W. E. Spicer, Physical Review B 5, 3017 (1972).
  * Leviton and Frey [2006] D. B. Leviton and B. J. Frey, Optomechanical Technologies for Astronomy 6273, 62732K (2006).
  * Steinlechner  _et al._ [2015] J. Steinlechner, I. W. Martin, J. Hough, C. Krüger, S. Rowan, and R. Schnabel, Physical Review D 91, 042001 (2015).
  * Fukuhara  _et al._ [1997] M. Fukuhara, A. Sanpei, and K. Shibuki, Journal of Materials Science 32, 1207 (1997).
  * Middelmann  _et al._ [2015] T. Middelmann, A. Walkov, G. Bartl, and R. Schödel, Physical Review B 92, 174113 (2015).
  * Adhikari  _et al._ [2020] R. X. Adhikari _et al._ , Classical and Quantum Gravity 37, 165003 (2020).
  * Murray  _et al._ [2015] P. G. Murray, I. W. Martin, K. Craig, J. Hough, R. Robie, S. Rowan, M. R. Abernathy, T. Pershing, and S. Penn, Physical Review D 92, 062001 (2015).
  * Birney  _et al._ [2018] R. Birney, J. Steinlechner, Z. Tornasi, S. MacFoy, D. Vine, A. S. Bell, D. Gibson, J. Hough, S. Rowan, P. Sortais, S. Sproules, S. Tait, I. W. Martin, and S. Reid, Physical Review Letters 121, 191101 (2018).
  * Steinlechner and Martin [2021] J. Steinlechner and I. W. Martin, Physical Review D 103, 042001 (2021).