†Daiki Tanabe<EMAIL_ADDRESS>figure table

# High-precision temperature monitoring system for room-temperature equipment
in astrophysical observations

Daiki Tanabea,b,c  Masaya Hasegawac,d  Masashi Hazumia,c,d,e,f  Nobuhiko
Katayamae  Shuhei Kikuchig  Adrian Leeh,i  Haruki Nishinoj  Satoru Takakurak
aDepartment of Particle and Nuclear Physics  School of High Energy Accelerator
Science  The Graduate University for Advanced Studies (SOKENDAI)  Tsukuba
Ibaraki 305-0801  Japan;
bCenter for High Energy and High Field Physics  National Central University
Taoyuan 32001  Taiwan;
cInstitute of Particle and Nuclear Studies (IPNS)  High Energy Accelerator
Research Organization (KEK)  Tsukuba  Ibaraki 305-0801  Japan;
dInternational Center for Quantum-field Measurement Systems for Studies of the
Universe and Particles (QUP)  High Energy Accelerator Research Organization
(KEK)  Tsukuba  Ibaraki 305-0801  Japan;
eKavli Institute for the Physics and Mathematics of the Universe (Kavli IPMU
WPI)  UTIAS  The University of Tokyo  Kashiwa  Chiba 277-8583  Japan;
fJapan Aerospace Exploration Agency (JAXA)  Institute of Space and
Astronautical Science (ISAS)  Sagamihara  Kanagawa 252-5210  Japan;
gYokohama National University  Yokohama  Kanagawa 240-8501  Japan;
hDepartment of Physics  University of California  Berkeley  CA 94720  USA;
iPhysics Division  Lawrence Berkeley National Laboratory  Berkeley  CA 94720
USA;
jResearch Center for the Early Universe  School of Science  The University of
Tokyo  Tokyo  Japan;111 Present affiliation: Japan Synchrotron Radiation
Research Institute (JASRI), 1-1-1, Kouto, Sayo-cho, Sayo-gun, Hyogo 679-5198,
Japan
kCASA & Department of Astrophysical & Planetary Sciences  University of
Colorado Boulder  Boulder  CO  80309  USA

###### Abstract

We present a precise thermometry system to monitor room-temperature components
of a telescope for radio-astronomy such as cosmic microwave background (CMB)
observation. The system realizes precision of 1 mK${\rm\sqrt{s}}$ on a
timescale of 20 seconds at 300 K. We achieved this high precision by tracking
only relative fluctuation and combining thermistors with a low-noise
measurement device. We show the required precision of temperature monitors for
CMB observation and introduce the performance of our thermometry system. This
precise room-temperature monitoring system enables us to reduce the low-
frequency noise in a wide range of radio-astronomical detector signals
observation and to operate a large detector array perform at its designed high
sensitivity.

###### keywords:

temperature monitor, thermistor, readout, mirror, cosmic microwave background

##  1 Introduction

Large detector arrays are now commonly used in astronomical and cosmological
observations. For millimeter and sub-millimeter imaging, the transition-edge
sensor (TES) bolometer allows high-precision measurements of radiation power.
The sensitivity of a single TES bolometer approaches the photon noise limit.
The multi-pixel TES bolometer array is thus becoming the standard for
improving sensitivity and expanding the field of view.

One of the major applications of the TES bolometer array is measuring the
cosmic microwave background (CMB). The odd-parity polarization pattern of the
CMB, which is called the B-mode, provides us with a probe of the inflation
hypothesis in cosmology on large angular scales [1]. Testing inflation with
the B-mode is, however, a challenge owing to the faint amplitude of the B-mode
on the degree scale. It is lower than the temperature anisotropy by five
orders of magnitude.

A standard strategy for measuring the B-mode precisely is improving
statistical sensitivity by increasing the number of detectors. Several ground-
based CMB experiments use $O(10,000)$ TES bolometers with designed array
sensitivities reaching $O(1)\ {\rm\mu K\sqrt{s}}$ [2, 3, 4, 5]. For example,
each receiver of the Simons Array (the POLARBEAR-2 telescope system) has 7,588
TES bolometers. Furthermore, the number of detectors for next-generation CMB
experiments is projected to be $O(500,000)$ in the late 2020s [6].

Despite the low noise of a single detector and the large number of detectors
in an array, noise correlated along multiple detectors cannot be mitigated by
averaging signals over the array. One of the sources of correlated noise is
the temperature fluctuation of the instruments. For example, the temperature
fluctuation of the focal plane, where the detectors are placed, modulates the
detector responsivity.

For a large array with more than $O(10,000)$ detectors, the temperature
fluctuation of instruments outside the cryostat is also not negligible
compared with the statistical sensitivity. There have been studies on the
pointing error caused by the structure deformation of a mirror and a telescope
due to the temperature variation [7, 8, 9]. The temperature fluctuation of the
mirror also modulates the polarizations of the reflection and the emission
from the mirror surface. This effect is seen by all detectors optically
coupled with the given area of the mirror. Another example is the temperature
fluctuation of readout electronics which changes the circuit gain. QUIET
experiment deployed temperature regulation for readout electronics [10]. In
recent experiments with large detector arrays, the fluctuation in readout
performance is correlated among detectors owing to multiplexing. These
correlated noises cannot be suppressed by averaging over detectors.

Whether we regulate temperature or correct the fluctuation effect in the
offline analysis, we have to monitor the temperature with sufficient
precision. If we observe a large angular scale with a relatively long
duration, it is necessary to track temperature fluctuation during
observations. There has been no report about room-temperature thermometry for
noise reduction from the third-generation CMB experiments. Several radio-
astronomical experiments have measured the room-temperature area for a
telescope foundation or mirror deformation, but the precision was insufficient
for reducing the noise of the CMB experiment [11, 12, 13].

In this paper, we show the requirements for precision of room-temperature
monitoring system and demonstrate that our thermometry system satisfies these
requirements. We focus on the mirror and the readout electronics as room-
temperature instruments on the basis of the design of the POLARBEAR-2
telescope system in the Simons Array. We assume the specifications of the
first telescope deployed in 2019 [14].

In section 2, we discuss the temperature coefficients between instrumental
temperature and a signal using the parameters of POLARBEAR-2, and then derive
the requirements for the temperature monitoring system. Section 3 introduces
the design concept and components of the temperature monitoring system. In
section 4, we summarize the monitor design and discuss the source of
systematic uncertainty in temperature measurement.

##  2 Requirements for the precise temperature monitoring

###  2.1 General formula for the temperature coefficient

To avoid degradation of the designed performance of a detector array, an
instrumental temperature fluctuation should be regulated or corrected so that
the excess noise stays below the statistical sensitivity of the detector
array. As we show later, the amplitude of the temperature fluctuation
corresponding to a statistical sensitivity of $O(10000)$-TES array is on the
order of $\mathrm{mK\sqrt{s}}$.

First we express the timestream of the detector signal as

$d(t)=s(t)+n(t)+\sum_{i}a_{i}T_{i}(t),$ (1)

where $s(t)$, $n(t)$, and $T_{i}(t)$ are functions of time indicating the
signal, noise, and temperature of the $i$-th instrument. Each instrumental
component has a temperature coefficient, $a_{i}$, which relates temperature to
the detector signal. It is a first-order perturbation expansion of the
temperature dependence. When the temperature fluctuation is correlated among
the detectors, we can treat $d(t)$ as an averaged timestream over an array.

The measured temperature $T^{\prime}_{i}(t)$ is

$T^{\prime}_{i}(t)=T_{i}(t)+\delta T_{i}(t),$ (2)

where $\delta T_{i}(t)$ is the uncertainty in temperature measurement. By
subtracting this measured temperature from the detector signal, we correct the
detector noise caused by temperature fluctuation. This operation can be
expressed via

$\begin{array}[]{lll}d^{\prime}(t)&\equiv&d(t)-\sum_{i}a_{i}T^{\prime}_{i}(t)\\\
&=&s(t)+n(t)-\sum_{i}a_{i}\delta T_{i}(t).\end{array}$ (3)

Thus, to avoid inducing excess noise, we require the uncertainty in
temperature measurement to be below the original array noise translated via
the temperature coefficient as

${\rm NSD}[\delta T_{i}]<\left|\frac{{\rm NSD}[n]}{a_{i}}\right|,$ (4)

where ${\rm NSD}[\delta T_{i}]$ and ${\rm NSD}[n]$ are the noise spectrum
densities of $\delta T_{i}(t)$ and $n(t)$, respectively. We focus on the
relative temperature fluctuation because the CMB observations do not involve
offset, but only measure the spatial variation of the CMB. Therefore, we
evaluate the detector array noise, uncertainty in thermometry, and actual
temperature fluctuation in terms of the noise spectrum density.

In the case of performing a temperature regulation, the $\delta T_{i}$ is
replaced with the temperature fluctuation after the regulation. In conclusion,
all instrumental temperature fluctuations have to be suppressed to stay below
the right side term of Eq. (4) by regulation or correction.

In the case of POLARBEAR-2, the designed array sensitivity of one receiver to
polarization is $4.1\ {\rm\mu K\sqrt{s}}$ [15]. In the following sections, we
use ${\rm NSD}[n]=1\ {\rm\mu K\sqrt{s}}$ as a fiducial array sensitivity to
polarization.

###  2.2 Room-temperature instruments

Here we assume a detection system that has bolometric detectors combined with
reflective optics, such as a CMB telescope. Many CMB experiments now use TES
bolometers for their small size, which allows a high-density array. Increasing
the number of detectors motivates the development of multiplexing techniques.
These techniques correlate readout noise with each other among detectors.
POLARBEAR-2 and SPT-3G have reflective optics and use the same frequency-
domain multiplexing (FDM) readout. See Appendix A for the bolometer principle
and FDM readout. AdvancedACT also has a reflective system, but the readout
method is time-domain multiplexing (TDM), which causes differences in
subdominant terms. Application to BICEP3 requires us drastic modifications of
the models because BICEP3 uses refractive optics and TDM readout.

Current CMB experiments aim to precisely measure the CMB polarization pattern.
One method of extracting polarization is pair-differences. This takes a
difference in signals from two detectors that are sensitive to two orthogonal
linear polarizations. This method is simple in concept, but the difference in
responsivities of the two detectors induces systematic error. Another way is
modulating polarization with a rotating half-wave plate before the light
arrives at the detector. This can extract polarization from only one detector
and can improve sensitivity at low frequency, but requires special care to
optically characterize the half-wave plate and stability of the rotation
encoder [16, 17]. To simplify the signal model, we consider excess noise
induced in the polarization by instrumental temperature fluctuation in the
pair-difference case.

The requirements are set at a certain frequency of interest. The frequency of
interest is determined by the scan speed of the telescope and the angular
scale of the targeted B-mode in the case of CMB observation. The scan speed of
POLARBEAR is $0.4^{\circ}/{\rm s}$ [18] and the inflationary B-mode has its
maximum amplitude at an approximately 3.6∘ scale. Therefore the temperature
needs to be stable over an interval of 9 seconds. We impose all requirements
at double this duration, namely, 50 mHz (20 seconds).

####  2.2.1 Readout electronics

We consider the temperature coefficient of a readout electronics system based
on the system used in POLARBEAR-2 and SPT-3G [19, 20]. As shown in Fig. 1, the
system consists of IceBoard motherboards, mezzanine boards, and
superconducting quantum interference device (SQUID) controller boards (SQCB)
[21]. A mezzanine board, which is stacked on an IceBoard, communicates with a
field-programmable gate array (FPGA) and processes signals to/from the
detectors. An SQCB applies bias current to the cooled SQUIDs and amplifies the
signal from the SQUIDs. Here we focus on the mezzanine board because the gain
fluctuation of an amplifier chain in the SQCB is compensated by digital active
nulling (DAN) feedback [22].

---
Figure 1:  Simplified diagram of the readout chain used in POLARBEAR-2.

The mezzanine board applies a constant bias voltage to the TESs, reads signals
from the SQUIDs, and sends a DAN feedback signal to the TESs to stabilize the
SQUIDs at their operating point. These signals fluctuate owing to the
temperature dependence of the resistors and ICs in the circuit.

In the pair-difference case, the fluctuation in bias voltage is the major
source of the excess polarization noise from the readout electronics. This
couples with the non-uniformity in the bias voltages applied to the pair
detectors and is induced in the polarization signal. Then the temperature
coefficient between the electronics temperature fluctuation and polarization
noise, $a_{\rm elec}$, can be expressed as

$a_{\rm elec}=\delta_{\rm bias}\frac{I_{\rm bias}}{g_{0}}A_{\rm bias},$ (5)

where $\delta_{\rm bias}$ is the fraction of non-uniformity in the bias
voltage, $g_{0}$ is a conversion from the Rayleigh–Jeans temperature to the
detector current, $I_{\rm bias}$ is the bias power applied to each TES, and
$A_{\rm bias}$ is a fractional coefficient between the circuit temperature and
bias voltage. The $\delta_{\rm bias}$ is estimated to be $0.5$% in POLARBEAR
[23]. The bias power in the temperature unit, $I_{\rm bias}/{g_{0}}$, is near
20 K for the Chilean observation. Assuming that the total temperature
dependence of bias generation circuit is $500$ppm/K [24], we obtain the
required precision of the readout temperature monitor as

${\rm NSD\left[T_{elec}\right]}<\left|\frac{{\rm NSD}[n]}{\displaystyle a_{\rm
elec}}\right|\sim 2.0\ {\rm mK\sqrt{s}}.$ (6)

####  2.2.2 Mirror

The polarized emission and reflection from the mirror can be estimated using
the formalism of Strozzi & McDonald (2000) based on the Fresnel law [25].
Their temperature dependence is simply linear under the Rayleigh–Jeans
approximation. The dominant polarization components are the polarized emission
from the mirror itself and polarization generated from the sky radiation when
it is reflected on the mirror surface. The total power of these components is
expressed in the Rayleigh–Jeans temperature unit as [26]

$P_{\rm mirror,pol}^{\rm
RJ}\simeq\displaystyle\sqrt{4\pi\rho\nu\epsilon_{0}}\frac{\sin^{2}\theta}{\cos\theta}\left(T_{\rm
mirror}-T_{\rm sky}\right),$ (7)

where $\rho$ is the resistivity of the mirror surface, $\nu$ is the frequency
of light, $\epsilon_{0}$ is the permittivity of vacuum, and $\theta$ is the
nutation angle between the receiver and the mirror. The quantities $T_{\rm
mirror}$ and $T_{\rm sky}$ are the temperatures of the mirror and sky,
respectively. Here we assumed that the permittivity of the air is the same as
that of vacuum. Note that $\rho$ has a slight temperature dependence.

The temperature dependence of this polarization is derived via

$\begin{array}[]{lll}a_{\rm mirror}&\equiv&\displaystyle\frac{dP_{\rm
mirror,pol}^{\rm RJ}}{dT_{\rm mirror}}\\\
&=&\displaystyle\sqrt{4\pi\rho\nu\epsilon_{0}}\frac{\sin^{2}\theta}{\cos\theta}\left\\{1-\frac{1}{2}\frac{d\ln\rho}{dT_{\rm
mirror}}\left(T_{\rm mirror}-T_{\rm sky}\right)\right\\},\end{array}$ (8)

Here we substitute parameters based on POLARBEAR-2. We fix the mean
temperature of the sky and the mirror as $T_{\rm sky}=10$ K [27] and $T_{\rm
mirror}=273$ K. Assuming the mirror made of aluminum, we use $\rho=2.417\times
10^{-8}\ \Omega\cdot{\rm m}$ at 273 K [28] and $d\ln\rho/dT_{\rm mirror}\sim
0.004$. The observation frequencies of light are $\nu=85.2$ and $143.9$ GHz
[14]. Because we are concerned with the upper limits of the temperature
coefficients, we use $\nu=143.9$ GHz as the light frequency. The angle between
the mirror surface and the detector plane is $\theta=32.5^{\circ}$[17]. Then
we obtain $a_{\rm mirror}\sim 0.032$%.

When this temperature coefficient is substituted into Eq.(4), the required
precision of the mirror temperature monitor is

${\rm NSD\left[T_{mirror}\right]}<\left|\frac{{\rm NSD}[n]}{a_{\rm
mirror}}\right|\sim 3.1\ {\rm mK\sqrt{s}}.$ (9)

For a Chilean telescope, the polarization signal from the mirror is maximally
modulated at sunrise and sunset because the temperature drift reaches 20 K on
the mirror surface at these times.

The required precision of the temperature monitor is on the order of ${\rm
mK\sqrt{s}}$ for both the readout electronics and the mirror.

We can evaluate the impacts of these temperature fluctuations to the detector
array noise using the temperature coefficients and measured temperature noise
spectra. The impacts of the room-temperature components are comparable. They
are also comparable to those of cold components such as the focal plane. The
temperature coefficient of the focal plane is thousands of times larger while
the magnitude of temperature fluctuation is on the order of ${\rm\mu
K\sqrt{s}}$ [29].

##  3 Design of a room-temperature monitoring system

###  3.1 Design concept

We realized a system that can monitor temperature fluctuation at 300 K with
precision on the ${\rm mK\sqrt{s}}$ level. Instead of facing difficulty in
achieving an accuracy of $O(1)$ mK with a room-temperature thermometer, we
focused on only tracking relative fluctuation. This motivates us to choose a
sensitive sensor even if it has uncertainty in absolute temperature. We
combined a low-noise measurement device with thermistors to push down the
precision limit determined by measurement noise.

We inserted a scanner to sense multiple thermometers on a large mirror
simultaneously. In addition, we ensured that the sampling rate of the
temperature monitors was higher than 100 mHz because we are concerned with
temperature fluctuation on the 50 mHz scale to study inflation.

We designed two types of room-temperature monitoring systems because the one
used for readout electronics needs special consideration for electromagnetic
compatibility. The temperature monitoring system for readout electronics is
shielded by a Faraday cage so as not to induce excess noise in the readout
through its electromagnetic emission.

###  3.2 Measurement device

We constructed a temperature monitoring system with an LTC2983 temperature
measurement system [30]. The LTC2983 can multiplex up to 10 channels with an
optional DC2210 terminal board and is controlled by a Linduino One
microcontroller.

We controlled the whole monitoring system with a Raspberry Pi 3B+ single-board
computer. The sampling rate was 0.7 Hz. The nominal uncertainty in readout
voltage was ${\rm\sim 1\ \mu V_{rms}}$, corresponding to 0.01 $\Omega$ when
the excitation current is 100 ${\rm\mu A}$.

We evaluated the system noise using a metal film resistor with a resistance of
3 ${\rm k\Omega}$ at 298 K [31]. The nominal temperature coefficient of this
resistor was 10 ppm/K, which was 10 times lower than the temperature
coefficients of the mirror and the readout electronics. We monitored the
resistance of this reference resistor with the LTC2983 system for an hour, and
calculated the noise spectrum density (NSD) for every 1000 seconds and
averaged the NSDs. Because the internal resistance of the LTC2983 itself has
temperature dependence, we monitored the internal resistance using a shorted
channel during the measurement and subtracted it from the data.

As shown in Fig. 2, the measured noise level around 50mHz was 0.04
${\rm\Omega\sqrt{s}}$, which corresponds to 0.3 ${\rm mK\sqrt{s}}$ with the
sensitivity of our thermistor. This was consistent with the expected noise
level based on the nominal readout uncertainty.

---
Figure 2:  Resistance measurement noise with a temperature measurement system
(LTC2983+DC2210).

For the mirror, we adopted a Keysight 34980A Multifunction Switch/Measure Unit
instead of LTC2983. It can multiplex up to 40 channels of two-wire systems
connected via an optional 34925A FET multiplexer and a 34925T terminal block,
with a 0.2-Hz sampling rate [32]. This feature is suitable for monitoring the
temperature of a wide area of the mirror with multiple thermometers.

We also measured the internal measurement noise of 34980A. The measurement
noise at 50mHz was ${\rm 0.06\ \Omega\sqrt{s}}$, corresponding to ${\rm 0.5\
mK\sqrt{s}}$ for our thermistor, which satisfies our requirement for
monitoring the mirror temperature. This was also consistent with the expected
noise level calculated from the nominal measurement uncertainty of the device.

###  3.3 Thermometer

We need to choose a highly sensitive thermometer to track a small fluctuation.
The sensor sensitivity can be expressed as the ratio between the change in a
sensor observable and the temperature drift.

We chose a thermistor for the thermometer. This is a ceramic resistor with a
nearly linear responsivity at room temperature. The thermistor response is
expressed in terms of the B-coefficient, which is determined via

$R=R_{0}\exp\left\\{B\left(\frac{1}{T}-\frac{1}{T_{0}}\right)\right\\},$ (10)

using two temperatures $T$ and $T_{0}$ and the corresponding resistances $R$
and $R_{0}$.

We tested a PS302J2 thermistor manufactured by Littelfuse Inc. [33]. Its
nominal resistance is $R=3000\ \Omega$ at $T=298$ K. The width of each sensor
is 2.4 mm. The sensor is guaranteed to operate in the temperature range from
$-80\ ^{\circ}$C to $+150\ ^{\circ}$C, which indicates adaptability to ambient
instruments in the major part of a year in the South Pole and Chile.

We calibrated it every 5 ${\rm{}^{\circ}C}$ in the range from $-30\ ^{\circ}C$
to $+30\ {\rm{}^{\circ}C}$ and obtained $B=3739$ in this temperature region.
This corresponds to a temperature coefficient of $\sim-123\ \Omega/{\rm K}$ at
298 K. This temperature coefficient and the measurement noise of LTC2983 give
us an expected minimum noise level of ${\rm 0.3\ mK\sqrt{s}}$.

Fig. 3 shows a comparison of the noise level with a platinum resistance
temperature detector (RTD) [34], thermocouple [35], and 1-wire thermometer
[36]. The thermistor had the lowest noise level as expected, which motivated
us to choose the thermistor as a temperature sensor.

---
Figure 3:  Temperature measurement noise of thermistor, platinum RTD, T-type
thermocouple, and 1-wire temperature sensor.

The noise level can be further pushed down by improving the calibration of the
T–R curve of the thermistors and applying a wind protection. As we show in
Sec. 3.4, we finally achieved a noise level of $0.6\pm 0.2\ {\rm mK\sqrt{s}}$.

###  3.4 Modification for on-site deployment

We adopted an LTC2983 as a measurement device for the readout electronics and
a Keysight34980A for the mirror. The thermistor was chosen as a thermometer.

We placed LTC2983 and its peripheral electronics in an aluminum box to shield
the electromagnetic emission from themselves. Its schematic and photograph are
shown in Fig. 4. There are two LTC2983 systems in the box. The box chassis is
grounded through the power ground.

---
Figure 4:  Temperature monitoring system for readout electronics. (a)
Schematic. (b) Photograph. The measuring system is located in a metal box for
RF shielding.

Figure 5 shows a schematic of wind protection and a sunshade applied to each
thermistor. First, we bonded a thermistor to a small and thin copper sheet
with an area of 10 mm square and thickness of 1 mm with varnish. Then we
covered it with silicone sealant and a thin reflective aluminum film. These
applications reduce the temperature fluctuation caused by wind, mitigate the
absolute temperature error caused by sunlight, and increase durability in the
severe environment of astronomical observatory.

---
Figure 5:  Scheme of the thermistor covered with silicone wind protection and
aluminum film.

We measured the temperature fluctuation of this covered thermistor with the
34980A in the laboratory. A group of thermistors was located in a thermal
chamber where hot wind circulates to control temperature. Half of them were
covered with silicone sealant and the other half were not. All of them were
thermally anchored to the metal floor of the thermal chamber through the
copper sheet. We calculated NSDs from 850-second cuts of time-ordered
temperature data and obtained ${\rm 2.5\pm 0.4\ mK\sqrt{s}}$ at 50 mHz for the
non-sealant group and ${\rm 0.8\pm 0.2\ mK\sqrt{s}}$ for the sealant group.
The NSDs are shown in Fig. 6. Without active temperature control by the wind,
the noise of the sealant group was further reduced to ${\rm 0.6\pm 0.2\
mK\sqrt{s}}$. This satisfies both requirements for the temperature monitor of
the mirror and the readout electronics, with the pair-difference model.

---
Figure 6:  Thermal noise of thermistors with- and without- silicone wind
protection, measured with 34980A+34925A. Each error bar is the standard
deviation of an 850-second window. The blue line is the averaged spectrum of
the non-sealant group. The peak around 65 mHz is due to the warm wind blowing
in the thermal chamber. The orange line is the averaged spectrum of the
sealant group.

We installed our warm temperature monitoring systems at the POLARBEAR-2
observation site in February 2019. We attached 12 thermistors on the back of
each mirror and 2 on each readout crate, and are continuously monitoring the
temperature of these warm instruments.

We measured the noise of the mirror temperature monitor at the observation
site. As shown in Fig. 7(a), 8 thermistors were radially attached to the outer
circumference of the mirror and 4 were around the center. We measured the
temperature of the thermistors from 2 a.m. to 8 a.m. (sunrise) and calculated
NSDs from 1000-second cuts. Figure 7(b) shows averaged NSDs of 4 thermistors
in the upper half including the left one viewed from the front, 4 in the lower
half including the right one, and 4 around the center. All thermistors have
temperature fluctuation below 3 ${\rm mK\sqrt{s}}$, indicating that they can
monitor the temperature with sufficiently low measurement noise at the site.
The variation of the fluctuation on the mirror could be induced by different
conditions of being exposed to wind.

---
Figure 7:  (a) Position of the 12 thermistors on the back of the primary
mirror of POLARBEAR-2. The grey region indicates the actual mirror area. (b)
Temperature fluctuation of thermistors on the primary mirror at night. NSDs of
4 thermistors in the upper half, 4 in the lower half, and 4 around the center
are averaged. Error bars are the standard deviation of every 4 NSDs.

##  4 Discussion

We should consider the temperature distribution of the target. For example,
temperature gradient between the back and front of a mirror, and that on the
front or back surface of a mirror cause errors in measuring absolute
temperature and relative temperature fluctuation.

We suggest reconstructing the temperature gradient on a surface using multiple
thermometers and subtracting the temperature fluctuation of a beam-projected
area multiplied by the temperature coefficient for every detector. We can
compensate for the temperature difference between the back and front of the
mirror by modeling the thermal conduction in the material and convection
caused by the air.

Another solution for treating a front-to-back temperature gradient is using
thermometers and wires thinner than the wavelength of the target light, and
directly attaching them to the front surface of the mirror. This is achievable
for observing radiation with a wavelength of over a millimeter. Embedding
sensitive thermometers in a circuit board is also effective for correcting the
temperature gradient around the readout electronics.

Our system for the mirror is designed to read many thermometers to monitor the
whole area of a large mirror. We tested 13 thermistors in this study. Our
measurement device for the mirror, a Keysight 34980A, can have its
multiplexing index increased to 320 by adding 34925A multiplexers. The
sampling rate, however, decreases when the number of thermometers increases.

We estimated the required precision of the temperature monitor was to be on
the order of millikelvin, but this can be loosened by configuring observation.
For example, a rotating half-wave plate can reduce the dependence on
temperature fluctuation of the readout electronics because it extracts the
polarization without uncertainty in the relative responsivity of detector
pairs.

##  5 Summary

We developed a system that can measure the temperature of room-temperature
instruments of a radio-astronomical telescope with high precision. We made one
system for the mirror and the other for the readout electronics. Both of our
monitoring systems consist of a low-noise measurement device, a multiplexer,
and thermistors. All components were commercial products, and the total cost
was approximately 4,000 USD. We achieved a noise spectrum density below ${\rm
1\ mK\sqrt{s}}$ per channel on a 20-second scale in measurements at 300 K in
the laboratory.

Both systems use a thermistor as a thermometer. This is the key to achieving
high precision owing to its high sensitivity and low cost. Because the concept
is tracking relative temperature fluctuation rather than measuring the
absolute value of temperature, we do not precisely calibrate the absolute
temperature scale of the thermometers. This enabled us to adopt a thermistor
in spite of its large absolute inaccuracy. The configuration of the thermistor
attachment was designed to block wind and direct sunlight, and to maintain a
fast response.

Our temperature monitoring system is applicable to a ground-based experiment
with a detector array with a noise level of $O(1)\ {\rm\mu{\rm K}\sqrt{s}}$.
This system will be an essential tool for observing subtle and continuous
signals, exemplified by the CMB.

##  6 Acknowledgement

This work was supported by JSPS KAKENHI grant numbers JP14J01662, JP18H01240,
JP18H04362, JP18J02133, and the SOKENDAI Short-term Research Abroad & Long-
term Internship Program. This research used resources of the Central Computing
System, owned and operated by the Computing Research Center at KEK. We
appreciate Dr. Neil Goeckner-Wald, who provided the basic idea of low-
frequency noise from instrumental temperature. Dr. Matt Dobbs, Dr. Joshua
Montgomery, Dr. Jean-François Cliche, and Dr. Tijmen de Haan provided detailed
information on the temperature dependence of the readout circuit. Dr. Grant
Teply, Dr. Nils Halverson, Dr. Brad Benson and Dr. Kevin Crowley informed us
about temperature monitoring systems in third-generation CMB experiments.
Additionally, Dr. Akito Kusaka, Dr. Aritoki Suzuki, Dr. Yuji Chinone and Dr.
John Groh gave us much advice on installing the temperature measurement
systems. Dr. Andrew May gave us insight in identifying noise sources by their
spectra. Mark Kurban from Edanz Group[37] edited a draft of this paper.

##  Appendix A Bolometer system in POLARBEAR-2

The TES bolometer is a high-resolution energy sensor operated in the
intermediate state between normal conduction and superconductivity, so called
the “transition edge” state. It is kept at this transition edge by applying
bias power to balance cooling and optical power. When a photon comes into the
heat absorber of the detector, the photon energy raises the temperature and is
detected as a sharp decrease in the bias current. This indicates that
fluctuation in the applied bias power directly affects the detector signal.
The bias voltage amplitude is typically 2 $\mu$V for POLARBEAR-2, which is
determined to be comparable to the optical load when the resistance of a TES
is $\sim 1\ \Omega$.

Increasing the number of TES bolometers motivates development of multiplexing
techniques. This is required for reducing the number of readout lines which
make the wiring complicated and transmit thermal load into the cryostat. In
the POLARBEAR-2 case, every 40 TES bolometers are read out by one line using
digital frequency-domain multiplexing (DfMUX) [38]. This method assigns a
different resonant frequency to each TES bolometer by connecting it to an LC
circuit. Signals from TES bolometers are modulated by the resonant
frequencies, accumulated on a readout line, amplified by superconducting
quantum interference devices (SQUIDs) on a 4 K stage, and finally demodulated
by electronics at the room temperature (Fig. 1). The gain fluctuation of this
readout chain makes the signal drift.

## References

  * [1] Uroš Seljak and Matias Zaldarriaga.  Signature of gravity waves in the polarization of the microwave background.  Phys. Rev. Lett., 78:2054–2057, Mar 1997.
  * [2] A. Suzuki et al.  The POLARBEAR-2 and the Simons Array Experiment.  J. Low. Temp. Phys., 184(3-4):805–810, 2016.
  * [3] S. W. Henderson, R. Allison, J. Austermann, et al.  Advanced ACTPol Cryogenic Detector Arrays and Readout.  Journal of Low Temperature Physics, 184(3-4):772–779, Mar 2016\.
  * [4] Z. Ahmed, M. Amiri, S. J. Benton, et al.  BICEP3: a 95GHz refracting telescope for degree-scale CMB polarization.  Millimeter, Submillimeter, and Far-Infrared Detectors and Instrumentation for Astronomy VII, Aug 2014.
  * [5] C M Posada, P A R Ade, Z Ahmed, et al.  Fabrication of large dual-polarized multichroic TES bolometer arrays for CMB measurements with the SPT-3G camera.  Superconductor Science and Technology, 28(9):094002, aug 2015.
  * [6] Kevork Abazajian, Graeme Addison, Peter Adshead, and others.  CMB-S4 Decadal Survey APC White Paper, 2019.
  * [7] J. E. Carlstrom, P. A. R. Ade, K. A. Aird, et al.  The 10 Meter South Pole Telescope.  Publications of the Astronomical Society of the Pacific, 123(903):568–581, may 2011.
  * [8] J. W. Fowler, M. D. Niemack, S. R. Dicker, et al.  Optical design of the Atacama Cosmology Telescope and the Millimeter Bolometric Array Camera.  Applied Optics, 46(17):3444, May 2007.
  * [9] R. J. Thornton, P. A. R. Ade, S. Aiola, et al.  THE ATACAMA COSMOLOGY TELESCOPE: THE POLARIZATION-SENSITIVE ACTPol INSTRUMENT.  The Astrophysical Journal Supplement Series, 227(2):21, Dec 2016\.
  * [10] C. Bischoff, A. Brizius, I. Buder, et al.  The Q/U Imaging ExperimenT Instrument.  The Astrophysical Journal, 768(1):9, 2013.
  * [11] Simon J. E. Radford Laura A. Snyder and Mark A. Holdaway.  ALMA Memo #314; Underground Temperature Fluctuations and Water Drainage at Chajnantor, Viewed in 2020.  http://legacy.nrao.edu/alma/memos/html-memos/abstracts/abs314.html.
  * [12] James Lamb and J. Forster.  MMA Memo #100; Temperature Measurements on BIMA 6-m Antennas - Part I: Backing Structure.  Viewed in 2020.  http://legacy.nrao.edu/alma/memos/html-memos/abstracts/abs100.html.
  * [13] E. V. Emelianov.  Analysis of thermal conditions of the 6-m BTA telescope elements and the telescope dome space.  Astrophysical Bulletin, 70(3):362?370, Jul 2015.
  * [14] Daisuke Kaneko et al.  Deployment of Polarbear-2A.  J. Low Temp. Phys., 199(3-4):1137–1147, 2020.
  * [15] Y. Inoue et al.  POLARBEAR-2: an instrument for CMB polarization measurements.  Proc. SPIE Int. Soc. Opt. Eng., 9914:99141I, 2016.
  * [16] Satoru Takakura et al.  Performance of a continuously rotating half-wave plate on the POLARBEAR telescope.  JCAP, 1705(05):008, 2017.
  * [17] Satoru Takakura.  Characterization of a continuous polarization modulator using a half-wave plate for measurements of degree-scale cosmic microwave background with the POLARBEAR experiment.  PhD thesis, Osaka University, 2017.
  * [18] S. Adachi et al.  A Measurement of the Degree Scale CMB $B$-mode Angular Power Spectrum with POLARBEAR.  10 2019.
  * [19] K. Hattori, Y. Akiba, K. Arnold, et al.  Development of Readout Electronics for POLARBEAR-2 Cosmic Microwave Background Experiment.  Journal of Low Temperature Physics, 184(1-2):512–518, Jan 2016\.
  * [20] A. N. Bender, P. A. R. Ade, A. J. Anderson, et al.  Integrated performance of a frequency domain multiplexing readout in the SPT-3G receiver.  In Wayne S. Holland and Jonas Zmuidzinas, editors, Millimeter, Submillimeter, and Far-Infrared Detectors and Instrumentation for Astronomy VIII, volume 9914, pages 321 – 331. International Society for Optics and Photonics, SPIE, 2016.
  * [21] K. Bandura, A. N. Bender, J. F. Cliche, et al.  ICE: a scalable, low-cost FPGA-based telescope signal processing and networking system, 2016.
  * [22] Tijmen de Haan, Graeme Smecher, and Matt Dobbs.  Improved Performance of TES Bolometers using Digital Feedback.  Proc. SPIE Int. Soc. Opt. Eng., 8452:84520E, 2012.
  * [23] P. A. R. Ade et al.  A Measurement of the Cosmic Microwave Background $B$-Mode Polarization Power Spectrum at Sub-Degree Scales from 2 years of POLARBEAR Data.  2017\.
  * [24] Author and J. Montgomery.  personal communication, 2021.
  * [25] David J, Strozzi, and Kirk T. McDonald.  Polarization Dependence of Emissivity, 2000.
  * [26] D. Barkats, C. Bischoff, P. Farese, et al.  Cosmic Microwave Background Polarimetry Using Correlation Receivers with the PIQUE and CAPMAP Experiments.  The Astrophysical Journal Supplement Series, 159(1):1–26, Jul 2005\.
  * [27] M. Hasegawa, O. Tajima, Y. Chinone, et al.  Calibration System with Cryogenically-Cooled Loads for CMB Polarization Detectors.  Rev. Sci. Instrum., 81:054501, 2011.
  * [28] P. D. Desai, H. M. James, and C. Y. Ho.  Electrical Resistivity of Aluminum and Manganese.  Journal of Physical and Chemical Reference Data, 13(4):1131–1172, 1984.
  * [29] Daiki Tanabe.  Precise room temperature monitoring for accurate measurements of degree-scale cosmic microwave background anisotropy at POLARBEAR-2.  PhD thesis, The Graduate University for Advanced Studies, 2020.
  * [30] Analog Devices.  LTC2983 Datasheet and Product Info, Viewed in 2021.  https://www.analog.com/en/products/ltc2983.html.
  * [31] Vishay Dale.  PTF Metal Film Resistors, Axial, High Precision, High Stability, Viewed in 2022.  https://www.vishay.com/ppg?31019l.
  * [32] Keysight Technologies.  34980A Multifunction Switch / Measure Unit, Viewed in 2021.  https://www.keysight.com/en/pd-429828-pn-34980A/.
  * [33] Littelfuse.  PS302J2-PS Series - Interchangeable Thermistors Leaded Thermistors, Viewed in 2021.  https://www.littelfuse.co.jp/products/temperature-sensors/leaded-thermistors/interchangeable-thermistors/standard-precision-ps/ps302j2.
  * [34] TE Connectivity.  NB-PTCO-011 : MEAS RTD Platinum Thin Film Element, Viewed in 2021.  https://www.te.com/usa-en/product-NB-PTCO-011.html.
  * [35] RS Components.  RS PRO Type T Thermocouple, 0.3mm Diameter, Viewed in 2021.  https://jp.rs-online.com/web/p/thermocouples/6212209/.
  * [36] Maxim Integrated.  Comparison of the DS18B20 and DS18S20 1-Wire Digital Thermometers, Viewed in 2021.  https://www.maximintegrated.com/jp/design/technical-documents/app-notes/4/4377.html.
  * [37] Edanz.  English Editing, Viewed in 2022.  https://jp.edanz.com/ac.
  * [38] Matt Dobbs, Eric Bissonnette, and Helmuth Spieler.  Digital Frequency Domain Multiplexer for Millimeter-Wavelength Telescopes.  IEEE Transactions on Nuclear Science, 55(1):21–26, 2008.

###### List of Figures

  1. 1  Simplified diagram of the readout chain used in POLARBEAR-2.
  2. 2  Resistance measurement noise with a temperature measurement system (LTC2983+DC2210).
  3. 3  Temperature measurement noise of thermistor, platinum RTD, T-type thermocouple, and 1-wire temperature sensor.
  4. 4  Temperature monitoring system for readout electronics. (a) Schematic. (b) Photograph. The measuring system is located in a metal box for RF shielding.
  5. 5  Scheme of the thermistor covered with silicone wind protection and aluminum film.
  6. 6  Thermal noise of thermistors with- and without- silicone wind protection, measured with 34980A+34925A. Each error bar is the standard deviation of an 850-second window. The blue line is the averaged spectrum of the non-sealant group. The peak around 65 mHz is due to the warm wind blowing in the thermal chamber. The orange line is the averaged spectrum of the sealant group.
  7. 7  (a) Position of the 12 thermistors on the back of the primary mirror of POLARBEAR-2. The grey region indicates the actual mirror area. (b) Temperature fluctuation of thermistors on the primary mirror at night. NSDs of 4 thermistors in the upper half, 4 in the lower half, and 4 around the center are averaged. Error bars are the standard deviation of every 4 NSDs.


# Spectral fit residuals as an indicator to increase model complexity

Anshuman Acharya Max Planck Institut für Astrophysik,
85748 - Garching, Germany.  Vinay L. Kashyap Center for Astrophysics $|$
Harvard & Smithsonian
60 Garden St. MS-70, Cambridge MA 02138, USA

###### Abstract

Spectral fitting of X-ray data usually involves minimizing statistics like the
chi-square and the Cash statistic. Here we discuss their limitations and
introduce two measures based on the cumulative sum (CuSum) of model residuals
to evaluate whether model complexity could be increased: the percentage of
bins exceeding a nominal threshold in a CuSum array (pctCuSum), and the excess
area under the CuSum compared to the nominal (p${}_{\textit{area}}$). We
demonstrate their use with an application to a Chandra ACIS spectral fit.

technique, statistical — software — spectral fitting, X-ray — Chandra

##  1 Introduction

Spectral fitting of X-ray data has been usually done by minimizing a statistic
like $\chi^{2}$ or cstat (Cash, 1979; Avni, 1976). The value of the statistic
is useful to evaluate the goodness of the fit, which is also used as a
stopping rule in evaluating model complexity. Fits are deemed acceptable if
$\chi^{2}/\nu<1+3\cdot\sqrt{\frac{2}{\nu}}$ (where $\nu$ is the degrees of
freedom in the fit) or $\Delta$cstat$<3\cdot\sigma_{\tt cstat}$ (where
$\sigma^{2}_{\tt cstat}$ is the estimated variance for the nominal model as
given by Kaastra 2017). When the fitting statistics or F-test values drop
below the threshold defined by the null distributions, it is commonly
acknowledged that further increases in model complexity (or the number of free
parameters in the fit) are statistically unsupported. However, these measures
are global measures of fit, and leave contiguous ranges in data space with
correlated deviations in the residuals. Here we introduce summary statistics
that identify the presence of such structures in the residuals and describe
how to use them to decide whether the model complexity should be increased.

We use a Chandra/ACIS-S spectrum of the corona of an exoplanet hosting star
(HD 179949, ObsID 6122; Acharya et al., 2023) to illustrate the measures (see
Figure 1). The data and the spectral fits made with different emission models
are shown in the left column of the Figure. The models increase in complexity
from top to bottom (see Acharya et al. 2023 for details), with none of the
statistical measures (cstat, pctCuSum, p${}_{\textit{area}}$) acceptable for
the simplest model, cstat being acceptable for the model in the middle row,
and all measures acceptable for the model in the bottom row. Below, we define
pctCuSum and p${}_{\textit{area}}$, describe their rationale, and how to
compute them and their null distributions.

##  2 CuSum: Definition, Calibration, and Significance

Spectral fit residuals can be expressed as

$\\{r_{i}(\hat{\theta})=(m_{i}(\hat{\theta})-c_{i}),~{}~{}~{}i=1,\ldots,N\\}\,,$

with $c_{i}$ the observed counts in bin $i$ of $N$ and $m_{i}$ the predicted
model counts in the same bin, generated for an astrophysical model defined by
the best-fit parameters $\hat{\theta}$. The cumulative sum (CuSum) at the
$j^{\rm th}$ bin is

${\rm CuSum}(\hat{\theta})_{j}=\sum_{i=1}^{j}r_{i}(\hat{\theta})\,.$ (1)

While the ordering of indices can be reversed, CuSum preserves the sequential
order of the bins, thus incorporating extra information ignored in $\chi^{2}$
or cstat statistics. The null distribution needed to evaluate the quality of
the CuSum of a best-fit can be built in one of two ways: $-$ generate $K$ mock
datasets using the fake_pha() function in Sherpa, and fit using the same
model, such that for each mock data set $k$, we obtain corresponding best-fit
parameters $\theta^{k}$; or $-$ obtain $K$ post burn-in draws
$\\{\theta^{k}\,,~{}k=1,\ldots,K\\}$ as iterations during a Markov Chain Monte
Carlo (MCMC) fit using the get_draws() tool in Sherpa (based on BLoCXS; van
Dyk et al., 2001). Each $\theta^{k}$ yields a CuSum array $\\{{\rm
CuSum}(\theta^{k})_{j}\,,~{}j=1,\ldots,N\\}$, and the resulting sample

$\\{{\rm CuSum}(\theta^{k})_{j}\,,~{}k=1,\ldots,K,~{}j=1,\ldots,N\\}$ (2)

defines the null distribution.

This null distribution characterizes the strength of the deviations present in
${\rm CuSum}(\hat{\theta})_{j}$. We construct two statistics that probe
deviation at different scales, pctCuSum to detect biases in the model
continuum and p${}_{\textit{area}}$ to detect presence/absence of narrow
lines.

  1. 1.

pctCuSum: We compute the equal-tail 90% point-wise range
$[C^{05}_{J},C^{95}_{J}]$ of ${\rm CuSum}(\theta^{k})_{j=J}$ for each detector
channel $J$ in the specified passband. Then we evaluate the percentage of
bins, pctCuSum, where the observed CuSum exceeds the 90% bounds, i.e., if $n$
channels have ${\rm CuSum}(\hat{\theta})_{J}<C^{05}_{J}$ or ${\rm
CuSum}(\hat{\theta})_{J}>C^{95}_{J}$, then
pct${}_{CuSum}=100\cdot\frac{n}{N}$. If pct${}_{CuSum}{\gtrsim}10$%, the model
used is considered inadequate and we suggest that more complex models with
more free parameters should be explored. If, on the other hand,
pct${}_{CuSum}{\ll}10$% this is a sign of overfitting, and thus less complex
models are favored (see the middle column of Figure 1).

  2. 2.

p${}_{\textit{area}}$: For each of the $k=1,\ldots,K$ CuSum arrays in the
null, we calculate the total area that falls outside the 5-95% bounds across
all the channels,

${\rm area}^{k}=\sum_{J}Z^{+}_{J}({\rm
CuSum}({\theta^{k}})_{J}-C^{95}_{J})+Z^{-}_{J}(C^{05}_{J}-{\rm
CuSum}({\theta^{k}})_{J})\,,$ (3)

where the indicator indices $Z^{+}=Z^{-}=1$ if the corresponding terms in the
summation are $+$ve, and $0$ otherwise. We then compute a similar area measure
for the best-fit spectrum,

${\rm area}(\hat{\theta})=\sum_{J}Z^{+}_{J}({\rm
CuSum}(\hat{\theta})_{J}-C^{95}_{J})+Z^{-}_{J}(C^{05}_{J}-{\rm
CuSum}(\hat{\theta})_{J})\,.$ (4)

The set $\\{{\rm area}^{k}\,,~{}k=1,\ldots,K\\}$ now forms a null distribution
of excess areas to which the value of ${\rm area}(\hat{\theta})$ can be
compared to obtain a $p$-value. The statistic p${}_{\textit{area}}$ is the
probability that ${\rm area}(\hat{\theta})>{\rm area}^{k}$. We consider the
best-fit model an inadequate fit if p${}_{area}{\ll}0.05$, as an indication
that a large deviation is present in ${\rm CuSum}(\hat{\theta})_{j}$ (see the
right column of Figure 1). This process is similar to the posterior predictive
$p$-value calibration procedure described by Protassov et al. (2002).

##  3 Application

In Figure 1, we show the results of the residuals analysis for three best-fit
models for an exemplar dataset of HD 179949 (Acharya et al., 2023). The top
row represents a single temperature APEC model with variable metallicity (1m).
The $\Delta{\tt cstat}=+3.3\cdot\sigma_{\tt cstat}$, indicates a poor fit.
This is supported by the pct${}_{CuSum}=40.6\%$ and a p-value = 0.0. The
middle row represents a 2-temperature model with metallicities scaling
together (2m). The $\Delta{\tt cstat}=+2.1\cdot\sigma_{\tt cstat}$, indicates
an acceptable model. The argument for improvements is solidified by
pct${}_{CuSum}=34.2\%$ and p-value = 0.0. The bottom row represents a
2-temperature model with abundances grouped by First Ionization Potential and
varying in several groups (2v). This gives $\Delta{\tt
cstat}=+1.3\cdot\sigma_{\tt cstat}$, pct${}_{CuSum}=13.5\%$ and p-value =
0.09. This model is therefore accepted as the best representation of the data.

|  |
---|---|---
|  |
|  |
Figure 1: Demonstrating the improvements in the measures of fit quality as
model complexity increases from top to bottom (see text). Each column
highlights a different statistic: $\Delta{\tt cstat}$ on left, pctCuSum in the
middle, and p${}_{\textit{area}}$ on right. The simplest model in the top row
fails in all measures; the $\Delta$cstat value in the middle row is borderline
acceptable, but the CuSum based measures are not; and in the bottom row all
the measures indicate that the fit is adequate. Increasing model complexity
further risks overfitting.

pctCuSum is sensitive to broad differences such as would be present in a
biased fit of the continuum, and p${}_{\textit{area}}$ to small scale changes
such as appear when spectral lines or abundance differences are mis-estimated.
Our method guards against overfitting to statistical fluctuations in the data
by specifying criteria for acceptability that are tied to the null
distributions. We note that our arguments and method are general, and can be
extended to any non-linear fitting scenario to supplement fitting procedures.

##  4 Code Availability

The code as described here is publicly available on Zenodo (Acharya & Kashyap,
2023). See https://github.com/anshuman1998/csresid for updates.

We thank Doug Burke (CXC/SDS) for suggestions with the code development.

## References

  * Acharya & Kashyap (2023) Acharya, A., & Kashyap, V. 2023, csresid: Spectral fit residuals as an indicator to increase model complexity, Zenodo, doi: 10.5281/zenodo.10395978
  * Acharya et al. (2023) Acharya, A., Kashyap, V. L., Saar, S. H., Singh, K. P., & Cuntz, M. 2023, ApJ, 951, 152, doi: 10.3847/1538-4357/acd054
  * Avni (1976) Avni, Y. 1976, ApJ, 210, 642, doi: 10.1086/154870
  * Cash (1979) Cash, W. 1979, ApJ, 228, 939, doi: 10.1086/156922
  * Kaastra (2017) Kaastra, J. S. 2017, Astronomy & Astrophysics, 605, A51, doi: 10.1051/0004-6361/201629319
  * Protassov et al. (2002) Protassov, R., van Dyk, D. A., Connors, A., Kashyap, V. L., & Siemiginowska, A. 2002, ApJ, 571, 545, doi: 10.1086/339856
  * van Dyk et al. (2001) van Dyk, D. A., Connors, A., Kashyap, V. L., & Siemiginowska, A. 2001, The Astrophysical Journal, 548, 224–243, doi: 10.1086/318656


# A lower bound for the Bogomolny and Schmit constant for random monochromatic
plane waves

Maxime Ingremeau and Alejandro Rivera

##  1 Introduction

Let $f:\mathbb{R}^{2}\rightarrow\mathbb{R}$ be the stationary, isotropic
planar centered almost surely continuous Gaussian field111We refer the reader
to [Abr97] for the definition and properties of random fields. with covariance
$\mathbb{E}[f(x)f(y)]=J_{0}(|x-y|)$ where
$J_{0}(r)=\frac{1}{2\pi}\int_{0}^{2\pi}e^{ir\sin(\theta)}d\theta$ is the $0$th
Bessel function of the first kind. We call this field the random monochromatic
plane wave. For each $R>0$, let us denote by $B(0,R)$ the ball of centre $0$
and of radius $R$, and by $N(R,f)$ the number of connected components of
$\mathbb{R}^{2}\backslash f^{-1}(\\{0\\})$ included in $B(0,R)$.

In [NS09], [NS16], Nazarov and Sodin showed that the limit

$\nu_{BS}=4\pi\times\lim_{R\rightarrow\infty}\frac{\mathbb{E}N(R,f)}{\pi
R^{2}}$ (1)

exists, and is positive. However, their method does not give an explicit value
for the constant $\nu_{BS}$, sometimes called the _Bogomolny-Schmit constant_
(or also Nazarov-Sodin constant). An equivalent definition for $\nu_{BS}$ is
that it is the limit of the average number of nodal domains for a random
spherical harmonic divided by the degree of this spherical harmonic. The
factor $4\pi$ can be interpreted as the area of the unit sphere. Bogomolny and
Schmit gave in [BS02] a highly heuristical argument, based on percolation,
which yielded the value

$\nu_{BS}\simeq 0.0624.$

However, numerical simulations carried out by Nastasescu ([Nas]), Konrad
([Kon12]) and Beliaev and Kereta ([BK13]) showed that $\nu_{BS}\simeq 0.0589$,
contradicting Bogomolny and Schmit’s prediction by a few percents.
Experimental measurements of $\nu_{BS}$ were also realised by Kuhl, Höhmann,
Stöckmann and Gnutzmann ([KHSG07]), using Berry’s conjecture ([Ber77]) that
high frequency eigenmodes in a chaotic cavity should locally behave like the
monochromatic random wave $f$. From the mathematical point of view, though,
very little is known regarding the value of $\nu_{BS}$. The best rigorous
lower bound so far was the one given in [Nas], which is of the order of
$10^{-319}$. The aim of this note is obtain a much better lower bound by
elementary means:

######  Theorem 1.1.

$\nu_{BS}\geq 1.39\times 10^{-4}.$

This bound is much smaller than the expected value of $\nu_{BS}$. However, our
method does not take into account all nodal domains, but only those which are
included in circles of radius $3.8$ (the first minimum of the Bessel function
$J_{0}$). After visual inspection of computer simulations we expect that these
are not very common. Aside from this, in our use of Lemma 2.2, we ignore the
fact that small, isolated nodal domains should be somewhat rare. We hope that
our methods can be used to count more general nodal domains, and to obtain
sharper lower bounds on $\nu_{BS}$.

Before proving Theorem 1.1 let us recall two basic facts about random
monochromatic plane waves.

  * •

The function $f$ almost surely satisfies $\Delta f+f=0$ where $\Delta$ is the
Laplace operator. This follows by applying the Laplace operator with respect
to $x$ to the expression $\mathbb{E}[f(x)f(y)]=J_{0}(|x-y|)$ because $J_{0}$
is also an eigenfunction of the Laplace operator.

  * •

In polar coordinates, $f$ can be given the simple expression

$f(r,\theta)=X_{0}J_{0}(r)-\sqrt{2}\sum_{n\geq
1}J_{n}(r)\left(X_{n}\cos(n\theta)+Y_{n}\sin(n\theta)\right)\,.$ (2)

where for each $n\in\mathbb{Z}$ we denote by $J_{n}$ the $n$th Bessel function
of the first kind, that is
$J_{n}(x)=\frac{1}{2\pi}\int_{0}^{2\pi}e^{in\theta-x\sin(\theta)}d\theta$, and
where $(X_{k})_{k\geq 0}$ and $(Y_{k})_{k\geq 1}$ are two families of centred
Gaussian random variables with unit variance and independent as a whole. An
explanation for this fact can be found for instance in Section 4.2 of [NS10].

### Acknowledgement

The authors would like to thank A. Deleporte for suggesting the proof of Lemma
2.2, as well as M. Sodin, S. Muirhead and M. McAuley for pointing out mistakes
in previous versions of this note.

The first author was funded by the Labex IRMIA, and partially supported by the
Agence Nationale de la Recherche project GeRaSic (ANR-13-BS01-0007-01).
The second author was funded by the ERC grant Liko No 676999.

##  2 Proof of Theorem 1.1

We now prove Theorem 1.1. The idea is to consider points such that $f$ does
not vanish on a certain circle around this point. We start by producing a
lower bound for the probability that $0$ is such a point.

######  Lemma 2.1.

Fix $r>0$ and write
$\Psi:t\mapsto\int_{t}^{+\infty}e^{-t^{2}/2}\frac{dt}{\sqrt{2\pi}}$. Let
$\mathcal{E}_{r}$ be the event that the zero set of $f$ does not intersect the
circle of radius $r$. Then, for each $T\in\mathbb{R}$,

$\mathbb{P}\left[\mathcal{E}_{r}\right]\geq
2\Psi(T)-\sqrt{2}r\Psi\left(\frac{T}{\sqrt{1-J_{0}(r)^{2}}}\right)\,.$

###### Proof.

In this proof we use the expression (2) for $f$. Fix $r>0$ and for each
$\theta\in[0,2\pi]$, let
$u(\theta)=X_{0}J_{0}(r)-f((r\cos(\theta),r\sin(\theta))$. We first fix
$x_{0}>0$ and try to estimate the probability that $u(\theta)$ crosses the
level $x_{0}J_{0}(r)$ when $\theta$ varies on the unit circle. Also,
throughout our calculations, we will use the following two Bessel function
identities:

$J_{-n}(x)=(-1)^{n}J_{n}(x);\,\sum_{n\in\mathbb{Z}}J_{n}(x)^{2}=1;\,\sum_{n\in\mathbb{Z}}n^{2}J_{n}(x)^{2}=\frac{x^{2}}{2}\,.$
(3)

These identities follow from the following classical formula (cf. [Wat22,
Chapter 2])

$\sum_{n\in\mathbb{Z}}t^{n}J_{n}(x)=e^{\frac{x}{2}(t-t^{-1})}$

by setting $t=e^{i\theta}$ and applying Parseval’s formula.
Now, observe that, since $f$ is stationary, isotropic and $(f,\nabla f)$ is
non-degenerate, $u$ is a stationary Gaussian field on the unit circle and the
field $(u,u^{\prime})$ is non-degenerate so that $u$ crosses the level
$x_{0}J_{0}(r)$ almost surely an even number of times. In particular, by
Markov’s inequality,

$\mathbb{P}\left[\exists\theta;\,u(\theta)=x_{0}J_{0}(r)\right]\leq\frac{1}{2}\mathbb{E}\left[\textup{Card}\\{\theta\in[0,2\pi];\,u(\theta)=x_{0}J_{0}(r)\\}\right]\,.$

To compute the right-hand side of this inequality, we apply the Kac-Rice
formula (see Theorem 6.2 from [AW09]). Define $\alpha(r)>0$ by
$\alpha(r)^{2}=\mathop{\mathrm{Var}}(u(\theta))=2\sum_{n\geq 1}J_{n}(r)^{2}$.
Then, by the first and second identity in (3), $\alpha(r)^{2}=1-J_{0}(r)^{2}$.
By the Kac-Rice formula,

$\mathbb{E}\left[\textup{Card}\\{\theta\in[0,2\pi];\,u(\theta)=x_{0}J_{0}(r)\\}\right]=\int_{0}^{2\pi}\mathbb{E}\left[|u^{\prime}(\theta)|\,|\,u(\theta)=x_{0}J_{0}(r)\right]\frac{e^{-\frac{x_{0}^{2}J_{0}(r)^{2}}{2\alpha(r)^{2}}}}{\alpha(r)\sqrt{2\pi}}d\theta\,.$

The fact that $u$ is stationary, implies first that the integrand is
independent of $\theta$ and second, that $u^{\prime}(\theta)$ is independent
of $u(\theta)$. Observe that
$\mathop{\mathrm{Var}}\left(u^{\prime}(\theta)\right)=2\sum_{n\geq
1}n^{2}J_{n}(r)^{2}$. By the first and third identities in (3),
$\mathop{\mathrm{Var}}(u^{\prime}(\theta))=\frac{r^{2}}{2}$. Moreover, if
$\xi$ is a centred Gaussian of variance one,
$\mathbb{E}\left[|\xi|\right]=\sqrt{\frac{2}{\pi}}$. Therefore,

$\mathbb{E}\left[\textup{Card}\\{\theta\in[0,2\pi];\,u(\theta)=x_{0}J_{0}(r)\\}\right]=\frac{\sqrt{2}r}{\alpha(r)}\exp\left(-\frac{x_{0}^{2}J_{0}(r)^{2}}{2\alpha(r)^{2}}\right)\,.$

On the other hand, we have

$\mathbb{P}[\neg\mathcal{E}_{r};\
X_{0}=x_{0}]=\mathbb{P}\big{[}\exists\theta\in[0,2\pi];\,u(\theta)=x_{0}J_{0}(r)\big{]}\leq\mathbb{E}\left[\textup{Card}\\{\theta\in[0,2\pi];\,u(\theta)=x_{0}J_{0}(r)\\}\right]$

In particular, for each $T\in\mathbb{R}$, the probability that $f$ has a zero
on the circle of radius $r$ centred at $0$ satisfies

$\mathbb{P}\left[\mathcal{E}_{r}\right]\geq\mathbb{E}\left[\left(1-\frac{r}{\sqrt{2}\alpha(r)}\exp\left(-\frac{X_{0}^{2}J_{0}(r)^{2}}{2\alpha(r)^{2}}\right)\right)\mathds{1}_{[|X_{0}|\geq
T]}\right]\,.$ (4)

Let $\Psi(t)=\mathbb{P}[X_{0}\geq
t]=\int_{t}^{+\infty}e^{-t^{2}/2}\frac{t}{\sqrt{2\pi}}$. Then, the right-hand
side is

$2\Psi(T)-2\frac{r}{\sqrt{2}\alpha(r)}\mathbb{E}\left[\exp\left(-\frac{X_{0}^{2}J_{0}(r)^{2}}{2\alpha(r)^{2}}\right)\mathds{1}_{[X_{0}\geq
T]}\right]\,.$

But for each $a>0$, a simple calculation shows that
$\mathbb{E}\left[\exp(-aX_{0}^{2})\mathds{1}_{[X_{0}\geq
T]}\right]=\frac{1}{\sqrt{1+2a}}\Psi\left(\sqrt{1+2a}T\right)$ so

$\mathbb{P}\left[\mathcal{E}_{r}\right]\geq
2\Psi(T)-\frac{\sqrt{2}r}{\sqrt{\alpha(r)^{2}+J_{0}(r)^{2}}}\Psi\left(\sqrt{1+J_{0}(r)^{2}/\alpha(r)^{2}}T\right)$

Replacing $\alpha(r)$ by its expression, we get, for any $T>0$ and any $r>0$,

$\mathbb{P}\left[\mathcal{E}_{r}\right]\geq
2\Psi(T)-\sqrt{2}r\Psi\left(\frac{T}{\sqrt{1-J_{0}(r)^{2}}}\right)\,.$ (5)

∎

Let $G_{r}=G_{r}(f)\subset\mathbb{R}^{2}$ be the (random) set of points
$x\in\mathbb{R}^{2}$ for which $f$ does not vanish on the circle of radius $r$
centred at $x$. If $x,y\in\mathbb{R}^{2}$, we will write $x\sim_{r}y$ if
$x,y\in G_{r}$ and $x,y$ belong to the same connected component of
$\mathbb{R}^{2}\backslash f^{-1}(\\{0\\})$. The next step of the proof is to
show that the connected components of $G_{r}$ are not too large.

######  Lemma 2.2.

Let $r_{1}$ and $r_{2}$ denote the first and second zeros of $J_{0}$ and let
$r\in]r_{1},r_{2}[$. Then, for each $x\in G_{r}$, the connected component of
$x$ in $\mathbb{R}^{2}\backslash f^{-1}(\\{0\\})$ is included in $B(x,r)$. In
particular, if $x,y\in G_{r}$ are such that $x\sim_{r}y$, then $|x-y|\leq r$.

######  Remark 2.3.

The result of this lemma may be optimal, deterministically speaking. But, we
expect equivalence classes of $G_{r}$ of diameter close to $r$ to be very
rare. We probably lose a large factor in this step. However, this intuition
seems difficult to quantify.

###### Proof of Lemma 2.2.

In this proof, we use the fact that $f$ satisfies $\Delta f+f=0$ almost
surely. Without loss of generality, we may suppose that $f(x)>0$. We claim
that $f(z)<0$ for all $z$ such that $|x-z|=r$. We already know that $f(z)$ has
constant sign for all $z$ such that $|x-z|=r$. Consider the function

$g(z):=\frac{1}{2\pi}\int_{0}^{2\pi}f(R_{\theta}z)\mathrm{d}\theta,$

where $R_{\theta}$ is the rotation of centre $x$ and of angle $\theta$. The
function $g$ satisfies $(\Delta+1)g=0$ (because $f(R_{\theta}\cdot)$ does for
each $\theta\in[0,2\pi]$), and it is radially symmetric around $x$. Therefore,
we must have

$g(z)=\lambda J_{0}(|z-x|)$

for some $\lambda\neq 0$. Since $f(x)>0$, we must have $\lambda>0$, and hence,
$f(z)<0$ for all $z$ such that $|x-z|=r$. Therefore, the connected component
of $\mathbb{R}^{2}\backslash f^{-1}(\\{0\\})$ containing $x$ is included in
$B(x,r)$ and the Lemma follows. ∎

We now finish off the proof of Theorem 1.1 by estimating the expected size of
$G_{r}$ using Lemma 2.1. Let $0<r_{1}<r_{2}$ be the first two zeros of the
Bessel function $J_{0}$ and fix $r\in]r_{1},r_{2}[$. Then, by Lemma 2.2, for
each $x\in G_{r}$, the equivalence class of $x$ has diameter at most $r$. By
the isodiametric inequality (see paragraph 10 of [BF87]) its area is no
greater than $\frac{\pi}{4}r^{2}$. Also, two different equivalence classes are
included in different connected components of $\mathbb{R}^{2}\backslash
f^{-1}(\\{0\\})$. Finally, if $R>r$ and $x\in B(0,R-r)$, then the connected
component of $\mathbb{R}^{2}\setminus f^{-1}\left(\\{0\\}\right)$ containing
$x$ is included in $B(0,R)$. Thus, for each $R>r$,

$\textup{Vol}\left(G_{r}\cap
B(0,R-r)\right)\leq\sum_{c}\textup{Vol}(c)\leq\frac{\pi}{4}r^{2}\times N(R,f)$

where the sum runs over the equivalence classes of $G_{r}$ intersecting
$B(0,R-r)$. Taking expectations, by stationarity, we get

$\textup{Vol}\left(B(R-r)\right)\mathbb{P}\left[0\in
G_{r}\right]\leq\frac{\pi}{4}r^{2}\times\mathbb{E}\left[N(R,f)\right]\,.$

Dividing by $\textup{Vol}\left(B(0,R)\right)=\pi R^{2}$ and letting
$R\rightarrow+\infty$, we get $\mathbb{P}\left[0\in
G_{r}\right]\leq\frac{1}{16}r^{2}\nu_{BS}$ which yields the following lower
bound for the Bogomolny-Schmit constant:

$\nu_{BS}\geq\frac{16\mathbb{P}[0\in G_{r}]}{r^{2}}\,.$

By Lemma 2.1, we have, for each $T>0$ and each $r>r_{0}$,

$\nu_{BS}\geq\frac{32}{r^{2}}\left[\Psi(T)-\frac{r}{\sqrt{2}}\Psi\left(\frac{T}{\sqrt{1-J_{0}(r)^{2}}}\right)\right]\,.$

Taking $r=3.8$ (the first minimum of $J_{0}$) and $T=3.35$ (the smallest $T$
for which in (4), the function whose expectation we compute is always
positive), we get

$\frac{32}{r^{2}}\geq 2.216\,;\ \frac{T}{\sqrt{1-J_{0}(r)^{2}}}\geq 3.659\,;\
\frac{r}{\sqrt{2}}\leq 2.69\,.$

We therefore obtain the announced lower bound

$\nu_{BS}\geq 1.39\times 10^{-4}\,.$

## References

  * [Abr97] P. Abrahamsen.  A review of Gaussian random fields and correlation functions.  Norsk Regnesentral/Norwegian Computing Center, 1997.
  * [AW09] J.-M. Azaïs and M. Wschebor.  Level sets and extrema of random processes and fields.  John Wiley & Sons, Inc., Hoboken, NJ, 2009.
  * [Ber77] M.V. Berry.  Regular and irregular semiclassical wavefunctions.  Journal of Physics A: Mathematical and General, 10(12):2083, 1977\.
  * [BF87] T. Bonnesen and W. Fenchel.  Theory of convex bodies. Transl. from the German and ed. by L. Boron, C. Christenson and B. Smith, with the collab. of W. Fenchel.  Moscow, Idaho (USA): BCS Associates. IX, 172 p. (1987)., 1987.
  * [BK13] D Beliaev and Z Kereta.  On the bogomolny–schmit conjecture.  Journal of Physics A: Mathematical and Theoretical, 46(45):455003, 2013.
  * [BS02] E. Bogomolny and C. Schmit.  Percolation model for nodal domains of chaotic wave functions.  Physical Review Letters, 88(11):114102, 2002.
  * [KHSG07] U. Kuhl, R. Höhmann, H.-J. Stöckmann, and S. Gnutzmann.  Nodal domains in open microwave systems.  Physical Review E, 75(3):036204, 2007.
  * [Kon12] Kyle Konrad.  Asymptotic statistics of nodal domains of quantum chaotic billiards in the semiclassical limit.  Senior Thesis, Dartmouth College, 2012.
  * [Nas] M Nastasescu.  The number of ovals of a real plane curve, Senior Thesis, Princeton 2011.  PhD thesis, Thesis and Mathematica code available at: http://www. its. caltech. edu/mnastase/Senior_Thesis. html.
  * [NS09] Fedor Nazarov and Mikhail Sodin.  On the number of nodal domains of random spherical harmonics.  American Journal of Mathematics, 131(5):1337–1357, 2009.
  * [NS10] F. Nazarov and M. Sodin.  Random complex zeroes and random nodal lines.  In Proceedings of the International Congress of Mathematicians 2010 (ICM 2010) (In 4 Volumes) Vol. I: Plenary Lectures and Ceremonies Vols. II–IV: Invited Lectures, pages 1450–1484. World Scientific, 2010.
  * [NS16] F. Nazarov and M. Sodin.  Asymptotic laws for the spatial distribution and the number of connected components of zero sets of Gaussian random functions.  Zh. Mat. Fiz. Anal. Geom., 12(3):205–278, 2016.
  * [Wat22] G Watson.  A Treatise on the Theory of Bessel Functions.  Cambridge University Press, 1922.