Silicon emissivity as a function of temperature
Marcio Constancio Jr
Divis ̃ao de Astrof ́ısica, Instituto Nacional de Pesquisas Espaciais, S ̃ao Jos ́e dos Campos,
SP, 12227-010,Brasil
Rana X. Adhikari
LIGO Laboratory, California Institute of Technology, Pasadena, CA 91125, USA
Odylio D. Aguiar
Divis ̃ao de Astrof ́ısica, Instituto Nacional de Pesquisas Espaciais, S ̃ao Jos ́e dos Campos,
SP, 12227-010,Brasil
Koji Arai
LIGO Laboratory, California Institute of Technology, Pasadena, CA 91125, USA
Aaron Markowitz
LIGO Laboratory, California Institute of Technology, Pasadena, CA 91125, USA
Marcos A. Okada
Divis ̃ao de Astrof ́ısica, Instituto Nacional de Pesquisas Espaciais, S ̃ao Jos ́e dos Campos,
SP, 12227-010,Brasil
Chris C. Wipf
LIGO Laboratory, California Institute of Technology, Pasadena, CA 91125, USA
Abstract
In this paper we present the temperature-dependent emissivity of a silicon
sample, estimated from its cool-down curve in a constant low temperature
environment (
∼
82
K
). The emissivity value follow a linear dependency in
the 120-260 K temperature range. This result is of great interest to the LIGO
Voyager gravitational wave interferometer project since it would mean that
no extra high thermal emissivity coating on the test masses would be required
in order to cool them down to 123 K. The results presented here indicate that
Preprint submitted to International Journal of Heat and Mass Transfer
April 27, 2020
arXiv:2004.06010v2 [physics.ins-det] 24 Apr 2020
bulk silicon itself can have sufficient thermal emissivity in order to cool the
200 kg LIGO Voyager test masses only by radiation in a reasonable short
amount of time (less than a week). However, it is still not clear if the natural
emissivity of silicon will be sufficient to maintain the LIGO Voyager test
masses at the desired temperature (123 K) while removing power absorbed
by the test masses. With the present results, a black coating on the barrel
surface of the test masses would be necessary if power in excess of 6 W
is delivered. However, the agreement we found between the hemispherical
emissivity obtained by a theory of semi-transparent Silicon and the obtained
experimental results makes us believe that the LIGO Voyager test masses,
because of their dimensions, will have effective emissivities around 0.7, which
would be enough to remove about 8.6 W (7.5 W) for a shield at 60 K (80K).
This hypothesis may be confirmed in the near future with new measurements.
Keywords:
LIGO, Silicon, thermal emissivity, gravitational waves
1. Introduction
Since the inauguration of gravitational wave astronomy in 2016, several
detections have been announced by the LIGO-VIRGO Scientific Collabora-
tion (LVC)[1–3]. From now on, a new way to observe the Universe is open,
with expectations ranging from completely new discoveries to the multi-
messenger astronomy with electromagnetic counterpart and even neutrinos.
Currently, advanced LIGO[4] and advanced Virgo[5], which are the 2
nd
generation of the LIGO and Virgo detectors, have regularly operated search-
ing for these signals. LIGO-India[6] and KAGRA[7] are expected to join this
journey within the next few years.
While the 2
nd
generation detectors are in operation, new detectors have
being planned. These include the Einstein Telescope (ET) [8] , LIGO Voyager[9]
and Cosmic Explorer detectors [10]. Between the last two, LIGO Voyager
brings great technological challenges since it is a cryogenic update in LIGO’s
current facilities.
In the Voyager version, Silicon-made mirror and suspensions will operate
at cryogenic temperatures due to some excellent properties, such as low me-
chanical loss in bulk silicon[11] and other optical properties at wavelengths
of 1.5 - 2.5
μ
m [12–14]. Also, Silicon has a zero crossing in thermal expan-
sion coefficient around 123 K, which can suppress thermoelastic noise at this
2
temperature [15–17]. The thermal expansion coefficient goes from -0.339 at
100 K to 2.618 at 300 K [15].
However, to keep both mirror and suspensions at this temperature, about
10 W of power deposited by the interferometer laser light must be extracted
from the test mass. This can be done through radiation if the emitting surface
has sufficiently high thermal emissivity. High emissivity coating such as DLC
- Diamond-Like Carbon coating [18] and Acktar Black coating [19] have been
proposed to improve the barrel’s thermal emissivity. The emissivity of dia-
mond like carbon was found to follow the linear fitting
DLC
= 0
.
3(
T/
300
K
),
which has lower values than the ones obtained for our sample. And the
emissivity of Acktar Black coating is
>
0.98 (in the 3-10
μ
m range) and
>
0.93 (in the 3-30
μ
m range). However, we need to know its emissivity in the
3-100
μ
m range or as a function of temperature to really be able to compare
its performance with our experimental results. Surface oxide layers can also
increase the emissivity properties. Relevant references are [20–24]. However,
any additional coating will deliver extra amounts of thermal noise to the
mirror.
In this paper we measure the emissivity of a bulk Silicon sample and show
that it may have a temperature-dependent thermal emissivity sufficiently
high in order to avoid the use of any extra coating on the barrel.
2. Experiment: Sample and Setup
In order to measure the total temperature-dependent emissivity of Silicon
(
(
T
)), a bulk Silicon sample was radiatively cooled down to 123 K while
having its temperature monitored by very thin thermocouples thermometers
(0.0799 mm in diameter).
The specimen used in this experiment is made of undoped, magnetic
Czochralski[25] grown silicon with a resistivity (as quoted by the manufac-
turer) of 4360 Ω.cm. It is 70.00 mm
×
30.00 mm
×
10.35 mm rectangular
parallelepiped shaped and weights 50.53 g. As shown in figure 1, it has a
rough finishing surface. The sample belongs to LIGO’s Caltech group but
the experiment was performed at Instituto Nacional de Pesquisas Espaciais
(INPE),Sao Jose dos Campos,Sao Paulo state, Brazil.
The experiment was performed in a cryostat under high vacuum (
<
10
−
7
mbar). This means convection and air conduction could be neglected. The
cryostat and a simplified diagram are shown in figure 2.
3
Figure 1: Silicon sample used in this experiment. The sample has not a mirror finishing
surface. T-type thermocouples and the Ti-6Al-4V wires are shown.
Figure 2: Cryostat. A picture of the real one (left) and a diagram with details of the
experiment (right).The stars show the position where the thermocouples were installed.
4
The sample was hang from the top of the dewar by four 0.29 mm thick,
20 cm long Ti-6Al-4V wires. Although there was a large gradient along the
wires, the heat transfer by them was 100-300 times lower than the dominant
heat transfer by radiation in the 123-300 K temperature range. Three T-type
thermocouples with 0.0799 mm in diameter (AWG40) were used to monitor
the temperature, two of them were attached to the Si sample and the other
in the top of the dewar (their position are shown with stars in figure 2).
In order to create a sudden cryogenic environment, the dewar was rapidly
immersed in a liquid nitrogen bath (
LN
2
). The amount of liquid was enough
to cool the dewar down and to guarantee that it would be surrounded by
LN
2
during the whole experiment. Finally, nine infrared baffles were used
to decrease the amount of heat coming from the top of the cryostat, which
was kept at room temperature. In order to check for reproducibility, the
experiment was performed twice.
3. Results and discussions
3.1. Analysis of heat transfer processes involved
Heat can be transferred from the Si sample to the cold dewar walls by
the following processes: radiation, gas thermal conductivity, gas thermal
convection, and metal thermal conductivity (by the suspension wires and
thermocouple wires). The values of power at 300 K and 123 K for these var-
ious processes are, respectively, the following: Radiation: 1.4 W (at 300K),
3.4
×
10
−
2
W (at 123K); gas thermal conductivity: 4.1
×
10
−
5
W (at 300K),
8.6
×
10
−
6
W (at 123K); gas thermal convection: negligible; thermal con-
ductivity of the suspension wires: 3.0
×
10
−
3
W (at 300K), 3.3
×
10
−
4
W (at
123K); thermal conductivity of the thermocouple wires: 1.8
×
10
−
3
W (at
300K), 4.1
×
10
−
4
W (at 123K). Therefore, radiation is by far the dominant
thermal transfer process. It is about 300 times and 46 times larger than all
the other thermal transfer processes at 300 K and 123 K, respectively. So in
the worse case (123K), all the other processes represent only about 2 percent.
3.2. Experimental results
As mentioned above, the experiment was performed twice in order to
check for reproducibility. The top of the dewar achieved 82.5(84.2) K in
about 6(4) minutes after the full immersion. This ensures that the whole
dewar was cold at this time. The final cooling down curves for runs 1(2) are
shown in figure 3. The inset shows the time when the sample crossed the
5
target 123 K temperature and shows the “noise” created by a digital reading.
Because the temperature reading is digital, the reading keeps going back and
forth, instead of showing a monotonic behavior.
Figure 3: Cooling down curve from run# 1 (left) and run# 2 (right). The inset shows the
time when the Silicon sample crossed the target 123 K temperature and its level of noise.
Considering that all the energy leaving the specimen is purely radiative
(convection, gas and wire conduction are negligible), the transient heat trans-
fer (or energy balance) equation can be written as:
m
Si
C
P
(
T
)
dT
dt
=
A
Si
(
T
)
σ
(
T
4
−
T
4
sh
)
(1)
where
(
T
) is the temperature-dependent total emissivity,
C
P
(
T
) is the
temperature-dependent heat capacity of Silicon[26],
σ
is the Stefan-Boltzmann
constant (= 5.6697
×
10
−
8
W.m
−
2
.K
−
4
),
T
Sh
is the temperature of the shield
(dewar) and
m
Si
,
A
Si
and T are, respectively, the mass, area and temperature
of the specimen (which is supposed to be uniform, given the slow cooling rate
and the high thermal conductivity of silicon). For a semi-transparent sample,
which is the case, as discussed in section 3.3, the use of the sample surface
area only makes sense for the purpose of defining an effective emissivity.
The emissivity can be calculated from the inverse problem by[27]:
(
T
) =
m
Si
C
P
(
T
)
dT
dt
A
Si
σ
(
T
4
−
T
4
sh
)
(2)
The derivative
dT
dt
can be obtained from the data in figure 3, making the
total temperature-dependent emissivity calculation feasible. However, since
6
the data is noisy, because the temperature reading is digital (the reading
keeps going back and forward), calculating the derivative from raw data can
affect its accuracy. So, a few different approaches were used to understand
the data. They are described in the following subsections.
3.2.1. First data analysis approach: Savitzky-Golay (SG) filtering
In the first approach, we used a Savitzky-Golay (SG) filter, which is a
digital filter based on a simplified least-squares fit convolution [28]. SG uses
a low-order polynomial in a smoothing window to fit the data. In this paper
the window size is 201 and the polynomial degree is 2. Scipy[29] package
in Python has a SG-function which gives the derivative of smoothed data
directly. The curve
dT
dt
obtained from this method was used as input in
equation 2. The emissivity calculated from this method is shown in figure 4
for both runs. The red curve is a linear regression of the data and the mean
between the both runs gives
(
T
) = 2
.
42
×
10
−
3
T
+ 1
.
22
×
10
−
1
.
Figure 4: Calculated emissivity for runs #1 and #2 by mean of SG filter. A mean value
between these linear regression gives
(
T
) = 2
.
42
×
10
−
3
T +1
.
22
×
10
−
1
.
3.2.2. Second approach: averaging the digital readings
Another approach we used was a calculation where a simplified aver-
age value was calculated from a set of points of the original array. A new
smoothed array
̄
T
was created from the original data (
T
), which is denoted
as (
T
1
,T
2
,T
3
,...T
i
). Every element
̄
T
j
results from the average value of
T
ranging from
i
−
n
to
i
+
n
of the original data, as follow:
̄
T
j
=
1
(2
n
+ 1)
k
=(
i
+
n
)
∑
k
=(
i
−
n
)
T
k
(3)
7
From this analysis, one can expect that
̄
T
j
has length equal to
len
(
T
)
−
2
n
.
In this analysis, n = 150, which means that the first 301 points (i-n to i+n)
from the original data were used to calculate T(151), which is the first point of
this array or
̄
T
1
(T(2) to T(302) will produce
̄
T
2
, and so on). The “derivative”
was calculated from this smoothed data by the difference
̄
T
(
i
−
n
)
−
̄
T
(
i
+
n
)
t
(
i
+
n
)
−
t
(
i
−
n
)
. As
before, the calculated
dT
dt
was used as input in equation 2.
The emissivity calculated through this method is shown in figure 5 for
both runs. The red curve is a linear regression of the data and the mean
value between the both runs gives
Si
(
T
) = 2.36
×
10
−
3
T + 1.14
×
10
−
1
.
Figure 5: Calculated emissivity for runs #1 and #2 by mean of a simplified average
smooth. A mean value between the linear regression of both runs gives
(
T
) = 2
.
36
×
10
−
3
T +1
.
14
×
10
−
1
.
3.2.3. Third approach: Curve fitting
Rather than smoothing (by averaging), as described above, another ap-
proach was tried to better understand the data and the results. Raw data was
fitted to theoretical curves before any calculation was performed. Two differ-
ent functions were used, a second order exponential function (
T
(
t
) =
a.e
−
b.t
+
c.e
−
d.t
+
e
)
and a 10
th
order polynomial function (
T
(
t
) =
a.t
10
+
b.t
9
+
c.t
8
+
d.t
7
+
e.t
6
+
f.t
5
+
g.t
4
+
h.t
3
+
i.t
2
+
j.t
+
k
).
From the fit, the dT/dt term was easily calculated and then
(
T
) was
obtained from equation 2. Figure 6 shows the
(
T
) versus T for both fittings
and for both runs. On the top, the result of the second order exponential fit
is shown for both runs. On the bottom, the same result is shown for the 10
th
order polynomial. For the exponential fit, the linear regression between
(
T
)
8
and T gives
(
T
) = 2
.
47
×
10
−
3
T +1
.
10
×
10
−
1
, On the other hand, the 10
th
order polynomial fit results in
(
T
) = 2
.
45
×
10
−
3
T +1
.
16
×
10
−
1
.
Figure 6: Calculated emissivity for runs #1 and #2 after curve fitting. On the top,
emissivity were calculated from a second order exponential fit and on the bottom they
were obtained from a 10
th
order polynomial fit. Polynomial fit is valid from 120-260K.
3.2.4. Forth approach: numerical differentiation by finite difference approx-
imation
Finally, the emissivity was also calculated from raw data, without any
curve fitting or averaging in the temperature array. However, to overcome
fluctuations in the data, the derivative took into account one point in the
(
i
−
n
)
th
position before and one point in the (
i
+
n
)
th
position after the
evaluated point. It means that for a given temperature
T
i
, the “derivative”
was calculated as the difference
∆
T
∆
t
=
T
i
−
n
−
T
i
+
n
t
i
+
n
−
t
i
−
n
.
This approach has been applied to a set of different values of
n
for runs #1
and #2. For all cases, the linear regression only had a maximum difference
of 2% and 6.5% for angular and linear coefficients, respectively. However, for
low
n
, the curve becomes noisier since it is more sensitive to the fluctuations
9
of temperature. Figure 7 shows the graphs for low and high values of
n
for
both runs. A mean value calculated between both runs using
n
= 150 gives
(
T
) = 2
.
43
×
10
−
3
T +1
.
21
×
10
−
1
.
Figure 7: Calculated emissivity for run #1 (left) and #2 (right) from raw data, without
curve fitting nor smoothing (by averaging) in the temperature array. The derivative was
calculated by taking into account one point in the (
i
−
n
)
th
position before and one point
in the (
i
+
n
)
th
position after the evaluated point.
3.2.5. Summary of results
All the values of emissivities with their respective coefficients of determi-
nation (
R
2
) for the various data analysis approaches are in Table 1.
In a general sense, the results obtained through all the methods described
here are in good agreement with each other, where the greater difference
between the average and the extremes was 6.5% and 20% for the angular
and linear coefficients, respectively. Although their good agreement, one of
them seem to have a best fit, which is the 10
th
order polynomial fit. The
emissivity values of this curve fit (the average between the one for run #1
10
Data analysis approach
(
T
)
R
2
Savitzky-Golay (SG)
filtering
Run #1
2
.
32
×
10
−
3
T + 1
.
42
×
10
−
1
0.8832197
Run #2
2
.
52
×
10
−
3
T + 1
.
01
×
10
−
1
0.9065968
Averaging the digital
readings
Run #1
2
.
27
×
10
−
3
T + 1.31
×
10
−
1
0.9654281
Run #2
2
.
44
×
10
−
3
T + 9.62
×
10
−
2
0.9608459
Curve fitting
(exp. 2
nd
order)
Run #1
2
.
42
×
10
−
3
T + 1.20
×
10
−
1
0.9519379
Run #2
2
.
51
×
10
−
3
T + 9.97
×
10
−
2
0.9734956
Curve fitting
(pol. 10
th
order)
Run #1
2
.
37
×
10
−
3
T + 1.32
×
10
−
1
0.9516301
Run #2
2
.
52
×
10
−
3
T + 1.00
×
10
−
1
0.9886076
Finite difference
approximation
Run #1
2
.
35
×
10
−
3
T + 1.37
×
10
−
1
0.8590187
Run #2
2
.
51
×
10
−
3
T + 1.04
×
10
−
1
0.8661256
Table 1: Summary of results
and the one for run #2) were used to determine a theoretical cooldown
curve, which was compared with the experimental data and shown in figure
8. This theoretical cooldown curve was calculated using equation 1 and
assuming arbitrary small steps of time (2 seconds). The new temperature
was calculated from the former one after this time step. This calculated
temperature would become the initial temperature for the next calculation
step and so on, everything starting from the initial temperature.
Figure 8 shows the experimental data for both experiments (blue and
green curves) and the calculated curve considering radiation and
(
T
) =
2
.
45
×
10
−
3
T
+ 1
.
16
×
10
−
1
. As we can see, this calculated function
(
T
)
greatly agrees with the data.
3.3. Simulation results
We performed Solidworks simulation to predict how long the cooling down
could last. In a first attempt, a constant emissivity of 0.5 was used. The blue
dotted curve in figure 9 show this temperature drop with time. A calculated
curve using only radiation, emissivity of 0.5 and equation 1, satisfactorily
predicted the same result, as shown in dashed black curve.
The experimental data for run #1 and #2 was also plotted in figure 9
to make a comparison easier. As can be seen, the cooling down time of the
Si sample from room temperature to 123 K was, coincidentally, about the
same as it had a constant emissivity of 0.5. Also, the experimental curve for
both runs is steeper in the beginning compared to the simulated 0.5 emissivity
11