# Age-Energy Trade-off in Status Update System with Wake-up Control

Jiajie Huang, Jie Gong This research was funded by the National Natural
Science Foundation of China grant number 62171481, the National Key Research
and Development Program of China grant number 2019YFE0114000, the Special
Support Program of Guangdong grant number 2019TQ05X150, the Natural Science
Foundation of Guangdong Province grant number 2021A1515011124, and the Science
and Technology Program of Guangzhou under Grant 202201011577. Guangdong Key
Laboratory of Information Security Technology,
School of Computer Science and Engineering, Sun Yat-sen University, Guangzhou,
510006, China
Email<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

In the status update system, the freshness of information is drawing more and
more attention. To measure the freshness of the data, age-of-synchronization
(AoS) is introduced. Since many communication devices are energy-constrained,
how to reduce energy consumption while keeping the timely response of data
needs to be carefully addressed. In this paper, we study the age-energy trade-
off problem in a single-source single-server scenario. We assume the server
enters a low-power sleep state when idle to save energy and consider three
wake-up policies. We adopt the stochastic hybrid system (SHS) method to
analyze the average AoS and average energy consumption under different
policies. The age-energy trade-off relationship under different parameters is
illustrated by numerical results.

###### Index Terms:

Age-of-synchronization, sleep-wakeup policy, stochastic hybrid system.

##  I Introduction

In the status update system, it is essential to know the state changes of the
source in time. As a result, many indicators to measure the freshness of
information have been proposed. Conventionally, age-of-information (AoI)[1]
has been widely studied, which is defined as the time elapsed since the latest
successfully accepted update was generated. Recently, age-of-synchronization
(AoS) is proposed to track whether the data is synchronized, which measures
the time elapsed since latest information at the receiver becomes
desynchronized [2]. Compared with the AoI metric, AoS is more suitable for
scenarios where the monitored data is updated less frequently, such as
databases, web crawling systems and error alarm systems. For the differences
in application scenarios between AoI and AoS, many performance studies on AoI
are not very suitable for AoS. Therefore, it is meaningful to study how to
obtain good AoS performance.

At the same time, due to the deployment of more and more communication
equipment, energy consumption has become one of the critical issues for
information and communication technology (ICT). Besides improving the
freshness of information, reducing the cost of communication energy
consumption also arouses widespread concern. However, improving information
freshness and saving energy consumption usually conflict with each other.
Intuitively, to keep the data fresh, it is necessary to process it in time,
which means that the equipment consumes more energy. Therefore, there is an
age-energy trade-off problem, and studying this problem has important guiding
significance in the state update system with limited energy consumption.

In recent years, there have been several studies on the trade-off between
information freshness and energy consumption. Energy harvesting sources were
considered in Ref. [3], while B. T. Bacinoglu et al. studied age optimal
strategies under infinite battery [4], unit battery [5], and finite battery
[6], respectively. The trade-off between energy consumption and AoI has
recently gained more attention [7, 8], and the age-energy trade-off has been
studied in error-prone channels [9] and in fading channels [10]. Ref. [11]
analyzed the age-energy trade-off in a state update system based on hybrid
automatic repeat request (HARQ). However, to the best of our knowledge, the
trade-off between AoS and energy consumption is an open problem and has not
been studied yet.

The SHS method is an effective and simple approach for analyzing data
freshness over networks, as proposed in [12]. It has been applied to various
scenarios and queuing models. For instance, [13] studies the average AoI of
each node in a single-source multi-hop status update system. The authors of
[14] consider a multi-source FCFS M/M/1 queuing model with an infinite queue
length. In [15], a multi-server and multi-source LCFS queuing model is
considered, with preemption in service adopted. The authors of [16] consider a
status update system with two sources and propose three packet management
strategies. The average AoI under each strategy is deduced by the SHS method.
However, energy consumption has not been taken into account in any of these
works. Therefore, this paper proposes to use the SHS method to analyze energy
consumption and AoS simultaneously.

In this paper, we focus on the trade-off between AoS and energy consumption.
We first introduce the system model and the definition of AoS. Secondly, a
sleep model is introduced and three wake-up policies are proposed. Then, the
SHS method is briefly introduced, and the average AoS and average energy
consumption of the three wake-up policies are analyzed by SHS. Finally,
numerical simulation shows the age-energy trade-off of different policies.

##  II System Model

In this work, we consider a single-source single-server real-time status
update system, as shown in Fig. 1. In this system, the status updates of the
source are generated randomly, which follows Poisson process with parameter
$\lambda$. The updates are sent to the server in the form of data packets. The
server processes the data packets at an exponentially distributed service rate
with parameter $\mu$, and sends them to the monitor. The preemption strategy
is adopted, that is, any newly generated packets from the source directly
preempt the ones being processed by the server.

Figure 1: System model. Figure 2: Sample AoS path. $\bullet$ indicates the
state update at the source and $\times$ indicates the state update at the
monitor is synchronized with the source.

We adopt the AoS [2] as the indicator of the freshness of the state
information of the physical process from the source. Fig. 2 describes a sample
path for source updates and monitor synchronization. Let’s define
${U_{1}},{U_{2}},\ldots{U_{k}}$ as the sequence of source update times,
${T_{1}},{T_{2}},\ldots{T_{k}}$ as the sequence of monitor refresh times. In
addition, we denote $N(t)$ as the number of source refreshes up to time $t$.
Formally, the definition of AoS is as follows.

###### Definition 1

_Let ${u}\left(t\right)$ denote the earliest time that the source gets a state
update since the last refresh of the monitor copy, i.e.,_

${u}\left(t\right)=\min\left\\{{\left.{{U_{k}}}\right|{U_{k}}>{T_{N\left(t\right)}}}\right\\}.$
(1)

_The AoS at time t is defined as_

$\Delta(t)=\max(t-u(t),0).$ (2)

_Note that if the monitor’s update is the same as the source, then
$\Delta(t)=0$._

It can be seen from (2) that whenever a new status update is generated from
the source, the monitor becomes unsynchronized with the source and the AoS of
the source starts to increase. The AoS value drops to zero when a new packet
is processed and received by the monitor, and remains at zero until a new
status update is generated. In our work, we adopt the average AoS of the
source over time, denoted by
$\bar{\Delta}={\lim_{t\to\infty}}E\left[{\Delta\left(t\right)}\right]$, as a
metric to evaluate the system performance.

Figure 3: Sleep model.

The energy consumption of the server is another important performance metric.
Therefore, we can adopt some sleep-wakeup policies to save energy when there
is no updates to be processed. The sleep model is illustrated in Fig. 3, where
the server has four states: _busy_ , _idle_ , _sleep_ and _wake-up_. In the
busy state, the server is processing data packets. After all the packets in
the server are processed, the server enters the idle state and operates at a
low power level. If any new packets arrive during the idle state, the server
immediately enters the busy state and starts processing them without delay or
energy cost. If no new packets arrive during the entire idle state, the server
enters the sleep state with extremely low power consumption, which can be
considered as zero. We assume that the duration of the idle state follows an
exponential distribution with mean $d$. The server transfers from the sleep
state to the wake-up state based on some wake-up policies. After a certain
amount of time in the wake-up state, the server turns to the busy state and
begins to process the packet. The time cost of the wake-up state is assumed to
follow an exponential distribution with mean $\theta$.

In our work, we consider the energy consumption rate ${\rm E}\left[P\right]$
as another performance metric. Specifically, we assume that the energy per
unit time consumed by the server in each state is denoted by
${P_{\textrm{B}}}$, ${P_{\textrm{I}}}$, ${P_{\textrm{S}}}$ and
${P_{\textrm{W}}}$ respectively. In general, ${P_{\textrm{B}}}$ is the
largest, ${P_{\textrm{S}}}$ is the smallest, and ${P_{\textrm{I}}}$ and
${P_{\textrm{W}}}$ are between them.

###  II-A Wake-up Policy

Different wake-up policies in queueing systems have been widely discussed in
literatures, such as N-policy [17], single-sleep [18], multiple-sleep [19]. In
this work, we consider these three wake-up policies, which are described in
detail below.

  * $\bullet$

N-policy: Under this wake-up policy, the server remains in the sleep state
until $N$ packets arrive. Once the ${N^{th}}$ packet arrives, it immediately
transfers to the wake-up state.

  * $\bullet$

Single-sleep: Under this wake-up policy, the server turns to the wake-up state
when it stays in the sleep state for a certain period of time. This period is
assumed to follow an exponential distribution with mean $s$. It is worth
noting that if there is no packet arrival during the wake-up state, the server
does not immediately turn into busy, but keep idle until a packet arrives.

  * $\bullet$

Multi-sleep: This wake-up policy can be seen as an extension of the Single-
sleep. In particular, the server firstly sleeps for a period of time after
entering the sleep state. If no data packets arrive during this period of
time, the server sleeps again for another period of time. The procedure
repeats until some data packets arrive. Then, the server transfers to the
wake-up state after the end of the current sleep period. Each time period is
also assumed to follow an exponential distribution with mean $s$.

In this paper, we aim to analyze the AoS and energy performance of the above
system with different wake-up policies. To study the problem with a unified
framework, we adopt the SHS method for analysis, which is introduced in the
next section.

##  III Analysis with Stochastic Hybrid System

In this section, we first briefly introduce the SHS method and show how to use
SHS to analyze the average AoS and the average energy consumption. Then, we
show the analytical results with three wake-up policies respectively.

###  III-A SHS method

####  III-A1 A Brief Introduction of SHS

SHS is a kind of stochastic dynamic system which combines continuous change
with discrete state variation, and the change of system structure is related
to some transformation rules[20]. The evolution of the discrete state is
determined by the transition or reset mapping, while the evolution of the
continuous state is determined by the stochastic differential equation. The
transition of discrete state is generally triggered by random events, and the
probability of transition at a given time depends on the continuous and
discrete components of the current SHS state. Therefore, SHS can be viewed as
a piecewise deterministic Markov process in a certain sense. According to the
definition of stochastic process, SHS can be expressed as

$\frac{{d{\mathbf{x}}\left(t\right)}}{{dt}}=f\left({q\left(t\right),{\mathbf{x}}\left(t\right),t}\right)+g\left({q\left(t\right),{\mathbf{x}}\left(t\right),t}\right)\frac{{d{\mathbf{z}}\left(t\right)}}{{dt}}.$
(3)

Where the discrete state is $q\left(t\right)\in\mathbb{Q}$, $\mathbb{Q}$ is a
discrete set. The continuous state is
${\mathbf{x}}\left(t\right)\in{\mathbb{R}^{n+1}}$.
${\mathbf{z}}\left(t\right)$ describes the process of independent Brownian
motion. Thus, there is a mapping
$f:\mathbb{Q}\times{\mathbb{R}^{n+1}}\times\left[{0,\infty}\right)\to{\mathbb{R}^{n+1}}$
and
$g:\mathbb{Q}\times{\mathbb{R}^{n+1}}\times\left[{0,\infty}\right)\to{\mathbb{R}^{\left({n+1}\right)\times
k}}$. And there is a set of transitions $L$, each $l\in L$ defines a discrete
transition/reset map
${\phi_{l}}:\mathbb{Q}\times{\mathbb{R}^{n+1}}\times\left[{0,\infty}\right)\to\mathbb{Q}\times{\mathbb{R}^{\left({n+1}\right)\times
k}}$. Therefore, the state transition is

$\left({q^{\prime}\left(t\right),{\mathbf{x^{\prime}}}\left(t\right)}\right)={\phi_{l}}\left({q\left(t\right),{\mathbf{x}}\left(t\right),t}\right).$
(4a) The corresponding transition intensity is
${\lambda^{\left(l\right)}}\left({q\left(t\right),{\mathbf{x}}\left(t\right),t}\right),\quad{\lambda^{\left(l\right)}}:\mathbb{Q}\times{\mathbb{R}^{n+1}}\times\left[{0,\infty}\right)\to\left[{0,\infty}\right).$
(4b)

When the system is in a discrete state, the continuous state evolves according
to (3). When the discrete state of the system changes from $q$ to
$q^{\prime}$, the continuous state jumps from $\mathbf{x}$ to
$\mathbf{x^{\prime}}$ according to (4a), and the frequency of the transition
is determined by (4b). In practice, the transition intensity is generally the
instantaneous rate at which the transition occurs.

####  III-A2 SHS for AoS and Energy Consumption

When using SHS to describe the AoS, the discrete state $q\left(t\right)$
represents the server occupancy, while the continuous state
$\mathbf{x}\left(t\right)$ represents the deterministic constant slope ramp
process. Therefore, for the general SHS model given in (3) and (4), we have

$f\left({q\left(t\right),{\mathbf{x}}\left(t\right),t}\right)={{\mathbf{b}}_{q}},$
(5a) $g\left({q\left(t\right),{\mathbf{x}}\left(t\right),t}\right)=0,$ (5b)
${\lambda^{\left(l\right)}}\left({q\left(t\right),{\mathbf{x}}\left(t\right),t}\right)={\lambda^{\left(l\right)}}{\delta_{{q_{l}},q}},$
(5c)
${\phi_{l}}\left({q\left(t\right),{\mathbf{x}}\left(t\right),t}\right)=\left({{{q^{\prime}}_{l}}\left(t\right),{\mathbf{x}}\left(t\right){{\mathbf{A}}_{l}}}\right).$
(5d)

In the Markov chain $q\left(t\right)$, each state $q\in\mathbb{Q}$ is a node
on the chain, and the transition between states $l$ is a directed edge
$({q_{l}},{q^{\prime}_{l}})$ with a transition rate of
${\lambda^{\left(l\right)}}$. The Cronecker function $\delta$ in (5c)
guarantees that the transition $l$ occurs only in the state $q_{l}$. When a
state transition occurs, the discrete state $q_{l}$ changes to the state
$q^{\prime}_{l}$, and the continuous state ${\mathbf{x}}\left(t\right)$ is
transformed according to the binary transfer reset mapping matrix
${{\mathbf{A}}_{l}}$:
${\mathbf{x^{\prime}}}\left(t\right)={\mathbf{x}}\left(t\right){{\mathbf{A}}_{l}}$.
In addition, according to (5a) and (5b), the evolvement of the continuous
state in each discrete state $q(t)=q$ is:

$\frac{{d{\mathbf{x}}\left(t\right)}}{{dt}}={{\mathbf{b}}_{q}}.$ (6)

Note that the evolution of AoS either increases with a slope of 1 or remains
unchanged at value 0. So in order to describe AoS, the continuous state vector
${\mathbf{x}}\left(t\right)$ in SHS requires only one dimension
${x_{0}}\left(t\right)$, corresponding ${{\mathbf{b}}_{q}}=b_{q}$, where
$b_{q}$ is a binary value. $b_{q}=1$ corresponds to the unit rate growth of
${x_{0}}\left(t\right)$ in discrete state $q$, and $b_{q}=0$ corresponds to
${x_{0}}\left(t\right)$ keeping constant in state $q$. Similarly, the values
of ${{\mathbf{A}}_{l}}={A}_{l}$ are only 0 and 1, ${A}_{l}=0$ means that the
update is completed and the AoS is reduced to 0, and ${A}_{l}=1$ means that
the AoS keeps increasing.

At the same time, when using SHS method to calculate the average energy
consumption and average AoS, it is necessary to calculate the stationary state
probability of Markov chain and the correlation vector between discrete state
${q}\left(t\right)$ and continuous state ${x_{0}}\left(t\right)$. Let
${{\mathbf{\pi}}_{q}}\left(t\right)$ denote the probability that the Markov
chain is in state $q$ , and
${{\mathbf{v}}_{q}}\left(t\right)={v_{{q_{0}}}}\left(t\right)$ denotes the
correlation between the discrete state and the continuous state. Therefore, we
can obtain

${{\mathbf{\pi}}_{q}}\left(t\right)=\Pr\left({q\left(t\right)=q}\right)={\rm
E}\left[{{\delta_{q,q\left(t\right)}}}\right].$ (7)
${{\mathbf{v}}_{q}}\left(t\right)={v_{{q_{0}}}}\left(t\right)={\rm
E}\left[{{x_{0}}\left(t\right){\delta_{q,q\left(t\right)}}}\right].$ (8)

Let $L_{q}$ represent the set of all transitions at state $q$, and
$L^{\prime}_{q}$ denote the set of transitions passed in at state $q$. One of
the basic assumptions for this kind of analysis is that Markov chain $q(t)$ is
ergodic. Under this assumption, the state probability vector
${\mathbf{\pi}}\left(t\right)=\left[{{\pi_{0}}\left(t\right)\ldots{\pi_{m}}\left(t\right)}\right]$
always converges to the only constant vector
${\mathbf{\bar{\pi}}}=\left[{{{\bar{\pi}}_{0}}\ldots{{\bar{\pi}}_{m}}}\right]$
which satisfies

${\bar{\pi}_{q}}\sum\nolimits_{l\in{L_{q}}}{{\lambda^{\left(l\right)}}}=\sum\nolimits_{l\in{L^{\prime}_{q}}}{{\lambda^{\left(l\right)}}}{\bar{\pi}_{{q_{l}}}},\quad
q\in\mathbb{Q},$ (9a) $\sum\nolimits_{q\in\mathbb{Q}}{{{\bar{\pi}}_{q}}=1},$
(9b)

In addition, if the Markov chain of the discrete state is ergodic and
stationary according to ${\mathbf{\bar{\pi}}}$, it has been shown in [12,
Theorem 4] that there must be a non-negative solution such that

${{\mathbf{\bar{v}}}_{q}}\sum\nolimits_{l\in{L_{q}}}{{\lambda^{\left(l\right)}}}={{\mathbf{b}}_{q}}{\bar{\pi}_{q}}+\sum\nolimits_{l\in{{L^{\prime}}_{q}}}{{\lambda^{\left(l\right)}}}{{\mathbf{\bar{v}}}_{{q_{l}}}}{{\mathbf{A}}_{l}},\quad
q\in\mathbb{Q},$ (10)

then the average AoS is given by

$\bar{\Delta}=\sum\nolimits_{q\in\mathbb{Q}}{{v_{{q_{0}}}}}.$ (11)

And the average energy consumption can be obtained by

${\rm
E}\left[P\right]=\sum\limits_{q\in\mathbb{Q}}{{{\bar{\pi}}_{q}}{P_{q}}},\quad{P_{q}}\in\left\\{{{P_{\textrm{B}}},{P_{\textrm{I}}},{P_{\textrm{S}}},{P_{\textrm{W}}}}\right\\}.$
(12)

###  III-B Analysis with SHS

In this part, we derive the average AoS and average energy consumption for
different wake-up policies with SHS.

####  III-B1 N-policy

In N-policy, the discrete state space of the Markov chain is
$\mathbb{Q}=\left\\{{B,ID,SL,1,2,\cdots,N}\right\\}$. In particular, $B$
refers to busy state, $ID$ refers to idle state, $SL$ refers to sleep state,
and $k\in\\{1,2,\cdots,N\\}$ refers to the state that a total of $k$ packets
have arrived during the sleep state. Note that $q\left(t\right)=N$
equivalently represents wake-up state, as the server immediately turns to this
state when the $N^{th}$ packet arrives. The continuous state degrades to a
scalar $x_{0}\left(t\right)=\Delta\left(t\right)$, which is the AoS of the
source.

Figure 4: State transition of N-policy.

The state transitions for the discrete state $q\left(t\right)$ is shown in
Fig. 4. The corresponding transitions of continuous state
${\mathbf{x}}\left(t\right)$ are summarized in the Table I. The transitions
are detailed as follows.

TABLE I: Table of Transitions for the Markov Chain in Fig. 4. $l$ | ${{q_{l}}\to{{q^{\prime}}_{l}}}$ | ${{\lambda^{\left(l\right)}}}$ | ${{\mathbf{x}}{{\mathbf{A}}_{l}}}$ | ${{{\mathbf{A}}_{l}}}$ | ${{{{\mathbf{\bar{v}}}}_{{q_{l}}}}{{\mathbf{A}}_{l}}}$
---|---|---|---|---|---
$1$ | ${B\to ID}$ | $\mu$ | ${\left[0\right]}$ | ${\left[0\right]}$ | ${\left[0\right]}$
$2$ | ${ID\to B}$ | $\lambda$ | ${\left[{{x_{0}}}\right]}$ | ${\left[1\right]}$ | ${\left[{{v_{ID0}}}\right]}$
$3$ | ${ID\to SL}$ | ${\frac{1}{d}}$ | ${\left[{{x_{0}}}\right]}$ | ${\left[1\right]}$ | ${\left[{{v_{ID0}}}\right]}$
$4$ | ${SL\to 1}$ | $\lambda$ | ${\left[{{x_{0}}}\right]}$ | ${\left[1\right]}$ | ${\left[{{v_{SL0}}}\right]}$
$5$ | ${1\to 2}$ | $\lambda$ | ${\left[{{x_{0}}}\right]}$ | ${\left[1\right]}$ | ${\left[{{v_{10}}}\right]}$
$\cdots$ | $\cdots$ | $\cdots$ | $\cdots$ | $\cdots$ | $\cdots$
$N+3$ | ${N-1\to N}$ | $\lambda$ | ${\left[{{x_{0}}}\right]}$ | ${\left[1\right]}$ | ${\left[{{v_{20}}}\right]}$
$N+4$ | ${N\to B}$ | ${\frac{1}{\theta}}$ | ${\left[{{x_{0}}}\right]}$ | ${\left[1\right]}$ | ${\left[{{v_{N0}}}\right]}$

  * -

$l=1$: When a packet completes service and is delivered to the monitor, the
server state changes from busy to idle. In this transition, the data on the
monitor side is synchronized with source. Thus, the AoS of source becomes
zero, i.e., $x^{\prime}_{0}=0$.

  * -

$l=2$: When a packet arrives in the idle state, the server turns to busy
state. In this case, the AoS of source remains the same, that is,
$x^{\prime}_{0}=x_{0}$, because the arrival of the packet does not reduce the
AoS until it is completely processed.

  * -

$l=3$: If no packets arrive during the idle state, the server turns to sleep
state. In this transition, $x^{\prime}_{0}=x_{0}$, because no packets are
processed.

  * -

$l=4,5,\cdots,N+3$: The state turns from state $l-4$ to $l-3$ when a new
packet arrives in the sleep state, where state 0 is equivalent to state $SL$.
This transition also does not change $x_{0}$.

  * -

$l=N+4$: When a total of $N$ packets arrive during sleep state, the server
immediately turns to wake-up state, and then turns to busy state after a
while. This transition also does not change $x_{0}$.

The evolution of ${\mathbf{x}}\left(t\right)$ is determined by the discrete
state $q\left(t\right)$. Specifically, when $q\left(t\right)=q$, we have

$\frac{{d{\mathbf{x}}\left(t\right)}}{{dt}}={{\mathbf{b}}_{q}}=\left\\{\begin{gathered}\left[1\right],q\in\\{1,2,\cdots,N,B\\}\hfill\\\
\left[0\right],q\in\\{ID,SL\\}\hfill\\\ \end{gathered}\right.$ (13)

The explanation for (13) is that $q\in\\{1,2,\cdots,N,B\\}$ means that there
are unprocessed packets in the system. Thus, the AoS grows at a unit rate.
$q\in\\{ID,SL\\}$ means that there are no packets in the system. Hence, the
source and the monitor are synchronized and the AoS remains constant at 0.

With the above conditions, we can calculate the average energy consumption and
average AoS by solving ${{{\bar{\pi}}_{q}}}$ and ${{\mathbf{\bar{v}}}_{q}}$.
Firstly, we use (9a) and (9b) to calculate the stationary probability vector
${\mathbf{\bar{\pi}}=\left[{\bar{\pi}_{B}},{\bar{\pi}_{ID}},{\bar{\pi}_{SL}},{\bar{\pi}_{1}},{\bar{\pi}_{2}},\cdots,{\bar{\pi}_{N}}\right]}$.
The matrix form of (9a) can be expressed as
${\mathbf{\bar{\pi}D}}={\mathbf{\bar{\pi}Q}}$, where $\mathbf{D}$ and
$\mathbf{Q}$ are given as

${\mathbf{D}}=diag\left[{\mu,\lambda+\frac{1}{d},\underbrace{\lambda,\lambda,\cdots,\lambda}_{Nelements},\frac{1}{\theta}}\right],$
${\mathbf{Q}}=\left[{\begin{array}[]{*{20}{c}}0&\mu&0&0&0&{}\hfil&0\\\
\lambda&0&{\frac{1}{d}}&0&0&{}\hfil&0\\\ 0&0&0&\lambda&0&\vdots&0\\\
0&0&0&0&\lambda&{}\hfil&0\\\
{}\hfil&{}\hfil&\cdots&{}\hfil&{}\hfil&\ddots&{}\hfil\\\
0&0&0&0&0&{}\hfil&\lambda\\\
{\frac{1}{\theta}}&0&0&0&0&\cdots&0\end{array}}\right].$

With ${\mathbf{\bar{\pi}D}}={\mathbf{\bar{\pi}Q}}$ and
$\sum\nolimits_{q\in\mathbb{Q}}{{{\bar{\pi}}_{q}}=1}$, we can obtain the
stationary probability of each state as

${{\bar{\pi}}_{B}}=A\cdot\frac{{1+d\lambda}}{\mu},$ (14a)
${{\bar{\pi}}_{ID}}=A\cdot d,$ (14b)
${{\bar{\pi}}_{SL}}={{\bar{\pi}}_{1}}=\cdots={{\bar{\pi}}_{N-1}}=A\cdot\frac{1}{\lambda},$
(14c) ${{\bar{\pi}}_{N}}=A\cdot\theta,$ (14d)

where
$A=\frac{1}{{\frac{N}{\lambda}+\frac{1}{\mu}+\theta+d\left({1+\frac{\lambda}{\mu}}\right)}}$.

By substituting (14) into (10), the value of ${\mathbf{\bar{v}}}_{q}$ is
obtained. Further substituting ${\mathbf{\bar{v}}}_{q}$ into (11), the average
AoS is given as

$\bar{\Delta}=\frac{{\frac{{1+d\lambda}}{\mu}+\frac{\theta}{\mu}+{\theta^{2}}+\frac{{N\left({N-1}\right)}}{{2{\lambda^{2}}}}+\frac{{N-1}}{\lambda}\left({\theta+\frac{1}{\mu}}\right)}}{{\frac{N}{\lambda}+\frac{1}{\mu}+\theta+d\left({1+\frac{\lambda}{\mu}}\right)}}.$
(15)

And the average energy consumption can be obtained by substituting (14) into
(12), which is expressed as

${\text{E}}\left[P\right]=\frac{{\frac{{1+d\lambda}}{\mu}{P_{\textrm{B}}}{+}d{P_{\textrm{I}}}{+}\frac{N}{\lambda}{P_{\textrm{S}}}{+}\theta{P_{\textrm{W}}}}}{{\frac{N}{\lambda}+\frac{1}{\mu}+\theta+d\left({1+\frac{\lambda}{\mu}}\right)}}.$
(16)

####  III-B2 Single-sleep policy

In single-sleep policy, the discrete state space of the Markov chain is
$\mathbb{Q}=\left\\{{SL,SL1,WK,WK1,B,ID0,ID}\right\\}$. To distinguish whether
the server is busy when a new packet arrives, we define two sets of states.
Specifically, the states $SL$, $WK$, $ID0$, $ID$ indicate that the server is
in sleep state, wake-up state, idle state after wake-up and idle state after
processing, respectively, without packets in the system. Note that the
difference between $ID0$ and $ID$ is that the server will turn to $SL$ if no
packets arrive during $ID$ state, but not in $ID0$. The states $SL1$, $WK1$,
$B$, indicate that the server is in sleep state, wake-up state, busy state,
respectively, all with packets in the queue or in processing. The continuous
state also degrades to a scalar $x_{0}\left(t\right)=\Delta\left(t\right)$.

Figure 5: State transition of single-sleep policy for a single source.

The state transition for the discrete state $q\left(t\right)$ is shown in Fig.
5. The corresponding transitions of continuous state
${\mathbf{x}}\left(t\right)$ are summarized in Table II. The transitions are
detailed as follows.

TABLE II: Table of Transitions for the Markov Chain in Fig.5. $l$ | ${{q_{l}}\to{{q^{\prime}}_{l}}}$ | ${{\lambda^{\left(l\right)}}}$ | ${{\mathbf{x}}{{\mathbf{A}}_{l}}}$ | ${{{\mathbf{A}}_{l}}}$ | ${{{{\mathbf{\bar{v}}}}_{{q_{l}}}}{{\mathbf{A}}_{l}}}$
---|---|---|---|---|---
$1$ | ${SL\to SL1}$ | $\lambda$ | ${\left[{{x_{0}}}\right]}$ | ${\left[1\right]}$ | ${\left[{{v_{SL0}}}\right]}$
$2$ | ${SL\to WK}$ | ${\frac{1}{s}}$ | ${\left[{{x_{0}}}\right]}$ | ${\left[1\right]}$ | ${\left[{{v_{SL0}}}\right]}$
$3$ | ${SL1\to WK1}$ | ${\frac{1}{s}}$ | ${\left[{{x_{0}}}\right]}$ | ${\left[1\right]}$ | ${\left[{{v_{SL10}}}\right]}$
$4$ | ${WK\to WK1}$ | $\lambda$ | ${\left[{{x_{0}}}\right]}$ | ${\left[1\right]}$ | ${\left[{{v_{WK0}}}\right]}$
$5$ | ${WK\to ID0}$ | ${\frac{1}{\theta}}$ | ${\left[{{x_{0}}}\right]}$ | ${\left[1\right]}$ | ${\left[{{v_{WK0}}}\right]}$
$6$ | ${WK1\to B}$ | ${\frac{1}{\theta}}$ | ${\left[{{x_{0}}}\right]}$ | ${\left[1\right]}$ | ${\left[{{v_{WK10}}}\right]}$
$7$ | ${ID0\to B}$ | $\lambda$ | ${\left[{{x_{0}}}\right]}$ | ${\left[1\right]}$ | ${\left[{{v_{ID0}}}\right]}$
$8$ | ${B\to ID}$ | $\mu$ | ${\left[0\right]}$ | ${\left[0\right]}$ | ${\left[0\right]}$
$9$ | ${ID\to B}$ | $\lambda$ | ${\left[{{x_{0}}}\right]}$ | ${\left[1\right]}$ | ${\left[{{v_{ID0}}}\right]}$
$10$ | ${ID\to SL}$ | ${\frac{1}{d}}$ | ${\left[{{x_{0}}}\right]}$ | ${\left[1\right]}$ | ${\left[{{v_{ID0}}}\right]}$

  * -

$l=1,4$: Packets arrive when the server is in the sleep or wake-up state. In
these transitions, the server turns into state SL1 or WK1. The AoS remains the
same, that is, $x^{\prime}_{0}=x_{0}$, because the arrival of the packet does
not reduce the AoS until it is completely processed.

  * -

$l=2,3$: After the server has been in sleep state for a period of time, it
turns into wake-up state. With this transition, the AoS of source does not
change, that is, $x^{\prime}_{0}=x_{0}$.

  * -

$l=5$: If no packets arrive during WK state, it turns into ID0 state to keep
idling and wait for a packet arrival. With this transition,
$x^{\prime}_{0}=x_{0}$.

  * -

$l=6$: After the server has been in WK1 state for a period of time, it turns
into busy state. With this transition, $x^{\prime}_{0}=x_{0}$.

  * -

$l=7$: Packets arrive when the server is in ID0 state. In this transition, the
server turns into busy state. With this transition, $x^{\prime}_{0}=x_{0}$.

  * -

$l=8$: A packet completes service, the server state changes from busy to idle.
In this transition, $x^{\prime}_{0}=0$.

  * -

$l=9,10$: The same as $l=2,3$ in Table I.

The evolution of ${\mathbf{x}}\left(t\right)$ is determined by the discrete
state $q\left(t\right)$. Specifically, when $q\left(t\right)=q$, we have

$\frac{{d{\mathbf{x}}\left(t\right)}}{{dt}}={{\mathbf{b}}_{q}}=\left\\{\begin{gathered}\left[1\right],q\in\\{SL1,WK1,B\\}\hfill\\\
\left[0\right],q\in\\{SL,WK,ID0,ID\\}\hfill\\\ \end{gathered}\right.$ (17)

The explanation for (17) is that $q\in\\{SL1,WK1,B\\}$ means that there are
unprocessed packets in the system. Thus, the AoS grows at a unit rate.
$q\in\\{SL,WK,ID0,ID\\}$ means that there are no packets in the system. Hence,
the source and the monitor are synchronized and the AoS at 0.

With the above conditions, we can calculate the average energy consumption and
average AoS by solving ${{{\bar{\pi}}_{q}}}$ and ${{\mathbf{\bar{v}}}_{q}}$.
Firstly, we use (9a) and (9b) to calculate the stationary probability vector
${\mathbf{\bar{\pi}}=\left[{\bar{\pi}_{SL}},{\bar{\pi}_{SL1}},{\bar{\pi}_{WK}},{\bar{\pi}_{WK1}},{\bar{\pi}_{B}},{\bar{\pi}_{B1}},{\bar{\pi}_{ID}}\right]}$.
The matrix form of (9a) can be expressed as
${\mathbf{\bar{\pi}D}}={\mathbf{\bar{\pi}Q}}$, where $\mathbf{D}$ and
$\mathbf{Q}$ are given as

${\mathbf{D}}=diag\left[{\lambda+\frac{1}{s},\frac{1}{s},\lambda+\frac{1}{\theta},\frac{1}{\theta},\lambda,\mu,\lambda+\frac{1}{d}}\right],$
${\mathbf{Q}}=\left[{\begin{array}[]{*{20}{c}}0&\lambda&\frac{1}{s}&0&0&0&0\\\
0&0&0&{\frac{1}{d}}&0&0&0\\\ 0&0&0&\lambda&{\frac{1}{\theta}}&0&0\\\
0&0&0&0&0&{\frac{1}{\theta}}&0\\\ 0&0&0&0&0&\lambda&0\\\ 0&0&0&0&0&0&\mu\\\
{\frac{1}{d}}&0&0&0&0&\lambda&0\end{array}}\right].$

With ${\mathbf{\bar{\pi}D}}={\mathbf{\bar{\pi}Q}}$ and
$\sum\nolimits_{q\in\mathbb{Q}}{{{\bar{\pi}}_{q}}=1}$, we can obtain the
stationary probability of each state as

${\bar{\pi}_{SL}}=\frac{1}{B}\left[{s\mu\lambda\left({\theta\lambda+1}\right)}\right],$
(18a)
${\bar{\pi}_{SL1}}=\frac{1}{B}\left[{{s^{2}}\mu{\lambda^{2}}\left({\theta\lambda+1}\right)}\right],$
(18b) ${\bar{\pi}_{WK}}=\frac{1}{B}\theta\mu\lambda,$ (18c)
${\bar{\pi}_{WK1}}=\frac{1}{B}\left[{\theta\mu{\lambda^{2}}\left({s+\theta+s\theta\lambda}\right)}\right],$
(18d) ${\bar{\pi}_{ID0}}=\frac{1}{B}\mu,$ (18e)
${\bar{\pi}_{B}}=\frac{1}{B}\left[{\lambda\left({d\lambda+1}\right)\left({s\lambda+1}\right)\left({\theta\lambda+1}\right)}\right],$
(18f)
${\bar{\pi}_{ID}}=\frac{1}{B}\left[{d\mu\lambda\left({s\lambda+1}\right)\left({\theta\lambda+1}\right)}\right],$
(18g)

where

$\begin{gathered}B=\mu+\lambda+ds\theta{\lambda^{4}}+\left({d+s+\theta}\right)\left({{\lambda^{2}}+\mu\lambda+s\theta\mu{\lambda^{3}}}\right)\hfill\\\
+\left({{s^{2}}+{\theta^{2}}+2s\theta+ds+d\theta}\right)\mu{\lambda^{2}}+\left({ds+s\theta+\theta
d}\right){\lambda^{3}}.\hfill\\\ \end{gathered}$

By substituting (18) into (10), the value of ${\mathbf{\bar{v}}}_{q}$ is
obtained. Further substituting ${\mathbf{\bar{v}}}_{q}$ into (11) and
$\mathbf{\bar{\pi}}$ into (12), the average AoS and average energy consumption
are given as

$\bar{\Delta}=\frac{{\lambda C}}{{\mu B}},$ (19)
${\text{E}}\left[P\right]=\frac{{s{P_{\textrm{S}}}{+}\theta{P_{\textrm{W}}}{+}({d{+}\frac{1}{{D}}}){P_{\textrm{I}}}{+}\frac{{d\lambda{+}1}}{\mu}{P_{\textrm{B}}}}}{{\frac{B}{\mu
D}}},$ (20)

where

$\begin{gathered}C=\left({{\mu^{2}}{\lambda^{2}}s\theta{+}\mu\lambda}\right)\left({{s^{2}}{+}s\theta{+}{\theta^{2}}}\right){+}{\mu^{2}}\lambda\left({{s^{3}}{+}{s^{2}}\theta{+}s{\theta^{2}}{+}{\theta^{3}}}\right)\hfill\\\
{+}\mu{\lambda^{2}}s\theta\left({s{+}\theta}\right){+}{\lambda^{3}}ds\theta{+}{\lambda^{2}}\left({ds{+}s\theta{+}\theta
d}\right){+}\lambda\left({d{+}s{+}\theta}\right){+}1,\end{gathered}$

$D=\lambda\left({s\lambda{+}1}\right)\left({\theta\lambda{+}1}\right)$.

####  III-B3 Multi-sleep policy

In Multi-sleep policy, the discrete state space of the Markov chain is
$\mathbb{Q}=\left\\{{SL,SL1,WK,B,ID}\right\\}$. In particular, all the states
have the same meaning as the ones with the same name in III-B2. The continuous
state also degrades to a scalar $x_{0}\left(t\right)=\Delta\left(t\right)$.

Figure 6: State transition of multi-sleep policy for a single sources.

The state transition for the discrete state $q\left(t\right)$ is shown in Fig.
6. The corresponding transitions of continuous state
${\mathbf{x}}\left(t\right)$ are summarized in Table III, which is very
similar to Table II. $l=1,3,4,5,6,7$ is the same as $l=1,3,5,8,9,10$ in Table
II. The main difference is the self-transition $l=2$, which means that the
server sleeps again after one sleep period if there is no packets arrival
during this period.

TABLE III: Table of Transitions for the Markov Chain in Fig.6. $l$ | ${{q_{l}}\to{{q^{\prime}}_{l}}}$ | ${{\lambda^{\left(l\right)}}}$ | ${{\mathbf{x}}{{\mathbf{A}}_{l}}}$ | ${{{\mathbf{A}}_{l}}}$ | ${{{{\mathbf{\bar{v}}}}_{{q_{l}}}}{{\mathbf{A}}_{l}}}$
---|---|---|---|---|---
$1$ | ${SL\to SL1}$ | $\lambda$ | ${\left[{{x_{0}}}\right]}$ | ${\left[1\right]}$ | ${\left[{{v_{SL0}}}\right]}$
$2$ | ${SL\to SL}$ | ${\frac{1}{s}}$ | ${\left[{{x_{0}}}\right]}$ | ${\left[1\right]}$ | ${\left[{{v_{SL0}}}\right]}$
$3$ | ${SL1\to WK}$ | ${\frac{1}{s}}$ | ${\left[{{x_{0}}}\right]}$ | ${\left[1\right]}$ | ${\left[{{v_{SL10}}}\right]}$
$4$ | ${WK\to B}$ | ${\frac{1}{\theta}}$ | ${\left[{{x_{0}}}\right]}$ | ${\left[1\right]}$ | ${\left[{{v_{WK0}}}\right]}$
$5$ | ${B\to ID}$ | $\mu$ | ${\left[0\right]}$ | ${\left[0\right]}$ | ${\left[0\right]}$
$6$ | ${ID\to B}$ | $\lambda$ | ${\left[{{x_{0}}}\right]}$ | ${\left[1\right]}$ | ${\left[{{v_{ID0}}}\right]}$
$7$ | ${ID\to SL}$ | ${\frac{1}{d}}$ | ${\left[{{x_{0}}}\right]}$ | ${\left[1\right]}$ | ${\left[{{v_{ID0}}}\right]}$

The evolution of ${\mathbf{x}}\left(t\right)$ is determined by the discrete
state $q\left(t\right)$. Specifically, when $q\left(t\right)=q$, we have

$\frac{{d{\mathbf{x}}\left(t\right)}}{{dt}}={{\mathbf{b}}_{q}}=\left\\{\begin{gathered}\left[1\right],q\in\\{SL1,WK,B\\}\hfill\\\
\left[0\right],q\in\\{SL,ID\\}\hfill\\\ \end{gathered}\right.$ (21)

The explanation for (21) is that $q\in\\{SL1,WK,B\\}$ means that there are
unprocessed packets in the system. Thus, the AoS grows at a unit rate.
$q\in\\{SL,ID\\}$ means there are no packets in the system. Hence, the source
and the monitor are synchronized and the AoS remains 0.

With the above conditions, we can calculate the average energy consumption and
average AoS by solving ${{{\bar{\pi}}_{q}}}$ and ${{\mathbf{\bar{v}}}_{q}}$.
Firstly, we use (9a) and (9b) to calculate the stationary probability vector
${\mathbf{\bar{\pi}}=\left[{\bar{\pi}_{SL}},{\bar{\pi}_{SL1}},{\bar{\pi}_{WK}},{\bar{\pi}_{B}},{\bar{\pi}_{ID}}\right]}$.
The matrix form of (9a) can be expressed as
${\mathbf{\bar{\pi}D}}={\mathbf{\bar{\pi}Q}}$, where $\mathbf{D}$ and
$\mathbf{Q}$ are given as

${\mathbf{D}}=diag\left[{\lambda+\frac{1}{s},\frac{1}{s},\frac{1}{\theta},\mu,\lambda+\frac{1}{d}}\right],$
${\mathbf{Q}}=\left[{\begin{array}[]{*{20}{c}}\frac{1}{s}&\lambda&0&0&0\\\
0&0&\frac{1}{s}&0&0\\\ 0&0&0&{\frac{1}{\theta}}&0\\\ 0&0&0&0&\mu\\\
{\frac{1}{d}}&0&0&\lambda&0\end{array}}\right].$

With ${\mathbf{\bar{\pi}D}}={\mathbf{\bar{\pi}Q}}$ and
$\sum\nolimits_{q\in\mathbb{Q}}{{{\bar{\pi}}_{q}}=1}$, we can obtain the
stationary probability of each state as

${{\bar{\pi}}_{SL}}=E\cdot\mu,$ (22a) ${{\bar{\pi}}_{SL1}}=E\cdot
s\mu\lambda,$ (22b) ${{\bar{\pi}}_{W}K}=E\cdot\theta\mu\lambda,$ (22c)
${{\bar{\pi}}_{B}}=E\cdot\lambda\left(1+d\lambda\right),$ (22d)
${{\bar{\pi}}_{ID}}=E\cdot d\mu\lambda,$ (22e)

where
$E=\frac{1}{{\mu+\lambda+d{\lambda^{2}}+d\mu\lambda+s\mu\lambda+\theta\mu\lambda}}$.

By substituting (22) into (10), the value of ${\mathbf{\bar{v}}}_{q}$ is
obtained. Further substituting ${\mathbf{\bar{v}}}_{q}$ into (11) and
$\mathbf{\bar{\pi}}$ into (12), the average AoS and average energy consumption
are given as

$\bar{\Delta}=\frac{{\lambda\left({{s^{2}}{\mu^{2}}{+}s\theta{\mu^{2}}{+}s\mu{+}{\theta^{2}}{\mu^{2}}{+}\theta\mu{+}d\lambda{+}1}\right)}}{{\mu\left({\mu{+}\lambda{+}d{\lambda^{2}}{+}d\mu\lambda{+}s\mu\lambda{+}\theta\mu\lambda}\right)}},$
(23)
${\text{E}}\left[P\right]=\frac{{\mu\left({s\lambda{+}1}\right){P_{\textrm{S}}}{+}\theta\mu\lambda{P_{\textrm{W}}}{+}\lambda\left({d\lambda{+}1}\right){P_{\textrm{B}}}{+}d\mu\lambda{P_{\textrm{I}}}}}{{\mu{+}\lambda{+}d{\lambda^{2}}{+}d\mu\lambda{+}s\mu\lambda{+}\theta\mu\lambda}}.$
(24)

##  IV Simulation Results

In this section, we present the trade-off between average AoS and average
energy consumption for the different policies through Monte Carlo simulations.
In the following, we assume that the energy consumed in each state in the
sleep model is ${P_{\textrm{B}}}=1$, ${P_{\textrm{I}}}={P_{\textrm{W}}}=0.5$,
${P_{\textrm{S}}}=0$, as per the assumptions made in our study.

Figure 7: Age-energy trade-off versus packets arrival rate
$\lambda\in[0.1,100]$, with $\mu=1$, $N=s=1$, $d=\theta=1$.

Fig. 7 depicts the trade-off between average AoS and average energy
consumption under different packet arrival rates $\lambda$. The simulation
results align well with the analytical ones, validating our theoretical
analysis. It is observed that the average energy consumption increases with
the increase of $\lambda$ in all the curves. This is due to the fact that the
higher the $\lambda$, the less likely the system is in a sleep state and the
more likely it is in a busy state, leading to larger energy consumption.
However, for average AoS, it does not always increase with the increase of
$\lambda$. The curves of single-sleep and multi-sleep exhibit a trend of
increasing first and then decreasing. This is because the increase of
$\lambda$ has a double effect on average AoS. When $\lambda$ is small, its
increase results in less time for information to be synchronized, thereby
increasing the average AoS. When $\lambda$ exceeds a certain threshold, its
increase increases the probability that the server transitions directly from
the idle state to the busy state, thereby accelerating the packet processing
process, leading to a reduction in the average AoS. Moreover, when $\lambda$
is large and tends to infinity, the results of the three policies are
indistinguishable. This is because the server does not enter sleep state but
always transitions directly from the idle state to the busy state.

Figure 8: The impact of $d$ on the age-energy trade-off with $\mu=1$, $N=s=1$,
$\theta=1$.

Then, we analyze the trade-off between average AoS and average energy
consumption under different system parameters in each wake-up policy. The
impact of idling time $d$ on the trade-off curve is shown in Fig. 8. As $d$
increases, the average energy consumption increases while the average AoS
decreases. This is because a higher $d$ leads to less time in the sleep state
and more time processing packets in the busy state.

Figure 9: The impact of $\theta$ on the age-energy trade-off with $\mu=1$,
$N=s=1$, $d=1$.

Fig. 9 shows the impact of wake-up time $\theta$ on the trade-off curve. The
average AoS increases as $\theta$ increases because it slows down the packet
processing, resulting in a larger average AoS. The average energy consumption,
however, shows different trends in different system loads. In the case of low
system load, where the arrival rate $\lambda$ is less than the service rate of
$\mu$, the server is in a low-power sleep state most of the time. Increasing
$\theta$ increases the time in the wake-up state, resulting in an increase in
average energy consumption. Conversely, in the case of heavy system load,
where the server is busy most of the time, increasing the sleep time reduces
the average energy consumption.

Figure 10: The impact of $N$ and $s$ on the age-energy trade-off with $\mu=1$,
$d=\theta=1$.

Fig. 10 illustrates the impact of sleep parameters ($N$ and $s$) on the age-
energy trade-off curve for different wake-up policies. The results show that
as $N$ or $s$ increases, the average energy consumption decreases, while the
average AoS increases in all policies. This is due to the direct effect of
increasing the sleeping time, which sacrifices the AoS to save energy.
Comparison of the three wake-up policies reveals that N-policy performs best
in terms of trade-off. Moreover, when the arrival rate of $\lambda$ is low, it
is not suitable to adopt N-policy with a higher $N$ value because it makes it
difficult to wake up when $N$ is large, leading to severe out-of-sync issues.
Furthermore, it is observed that as the arrival rate of $\lambda$ increases,
the trade-off curves of single-sleep policy and multi-sleep policy almost
coincide because fast packet arrival makes it unlikely that there will be no
packet arrival within a sleep duration $s$, and the effect of multi-sleep
policy is almost equivalent to that of single-sleep policy.

##  V Conclusion

In this paper, we investigate a status update system comprising a source and a
server, and propose a sleep model to conserve energy when the server is idle.
We introduce three wake-up policies to activate the server and derive explicit
expressions of the average AoS and average energy consumption using the SHS
method for each policy. We demonstrate the trade-offs between average AoS and
average energy consumption under different system parameters via Monte Carlo
simulation. Our simulation results show that, for a heavy system load, the
sleep model effectively limits the growth of AoS. For a light system load, a
small wake-up time $\theta$ can reduce both age and energy consumption. The
trade-off scope that N-policy provides widens as the system load decreases.
For a high system load, N-policy offers the best age performance, while
single-sleep and multi-sleep policies tend to be comparable. In future work,
we plan to explore the expansion of the system to multiple sources.

## References

  * [1] S. Kaul, R. Yates and M. Gruteser, “Real-time status: How often should one update?,” in Proc. IEEE Conf. Comput. Commun. (INFOCOM), Mar. 2012, pp. 2731–2735.
  * [2] J. Zhong, R. D. Yates and E. Soljanin, “Two Freshness Metrics for Local Cache Refresh,” in Proc. IEEE Int. Symp. Inf. Theory (ISIT), Jun. 2018, pp. 1924-1928.
  * [3] R. D. Yates, “Lazy is timely: Status updates by an energy harvesting source,” in Proc. IEEE Int. Symp. Inf. Theory (ISIT), Jun. 2015, pp. 3008–3012.
  * [4] B. T. Bacinoglu, E. T. Ceran, and E. Uysal-Biyikoglu, “Age of information under energy replenishment constraints,” in Proc. Inf. Theory Appl. Workshop (ITA), Feb. 2015, pp. 25–31.
  * [5] B. T. Bacinoglu and E. Uysal-Biyikoglu, “Scheduling status updates to minimize age of information with an energy harvesting sensor,” in Proc. IEEE Int. Symp. Inf. Theory (ISIT), Jun. 2017, pp. 1122–1126.
  * [6] B. T. Bacinoglu, Y. Sun, E. Uysal, and V. Mutlu, “Optimal status updating with a finite-battery energy harvesting source,” J. Commun. Netw., vol. 21, no. 3, pp. 280–294, Jun. 2019.
  * [7] M. Xie, J. Gong, X. Jia, X. Ma, “Age and Energy Tradeoff for Multicast Networks With Short Packet Transmissions,” IEEE Trans. Commun., vol. 69, no. 9, pp. 6106 - 6119, Sept. 2021.
  * [8] M. Xie, J. Gong, Q. Wang, X. Ma, “Age and Energy Analysis for LDPC Coded Status Update with and without ARQ,” IEEE Internet Things J., vol. 7, no. 10, pp.10388 - 10400, Oct. 2020.
  * [9] J. Gong, J. Zhu, X. Chen, X. Ma, “Sleep, Sense or Transmit: Energy-Age Tradeoff for Status Update with Two-Thresholds Optimal Policy,” IEEE Trans. Wireless Commun., vol. 21, no. 3, pp. 1751-1765, Mar. 2022.
  * [10] H. Huang, D. Qiao, and M. C. Gursoy, “Age-energy tradeoff in fading channels with packet-based transmissions,” in Proc. IEEE Conf. Comput. Commun. (INFOCOM) Workshops, Jul. 2020, pp. 1–6.
  * [11] J. Huang, J. Gong and W. Zhan, “Energy-Age Tradeoff in Status Update Communication Systems Based on HARQ,” in Proc. IEEE WCSP, Oct. 2021, pp. 1–5.
  * [12] R. D. Yates and S. K. Kaul, “The age of information: Real-time status updating by multiple sources,” IEEE Trans. Inf. Theory, vol. 65, no. 3, pp. 1807–1827, Mar. 2019.
  * [13] S. K. Kaul and R. D. Yates, “Age of information: Updates with priority,” in Proc. IEEE Int. Symp. Inf. Theory (ISIT), Jun. 2018, pp. 2644–2648.
  * [14] S. K. Kaul and R. D. Yates, “Timely updates by multiple sources: The M/M/1 queue revisited,” in Proc. 54th Annu. Conf. Inf. Sci. Syst. (CISS), Mar. 2020, pp. 1–6.
  * [15] A. Javani, M. Zorgui, and Z.Wang, “Age of information in multiple sensing,” in Proc. IEEE Global Commun. Conf. (GLOBECOM), Dec. 2019, pp. 9–13.
  * [16] M. Moltafet, M. Leinonen, and M. Codreanu, “Average AoI in multi-source systems with source-aware packet management,” IEEE Trans. Commun., vol. 69, no. 2, pp. 1121–1133, Feb. 2021.
  * [17] J. C. Ke, “The optimal control of an m/g/1 queueing system with server vacations, startup and breakdowns,” Comp & Ind. Eng., vol. 44, no. 4, pp. 567–579, 2003.
  * [18] W. Liu, Y. Ma, and J. Li, “Equilibrium threshold strategies in observable queueing systems under single vacation policy,” Appl. Math. Model., vol. 36, no. 12, pp. 6186–6202, 2012.
  * [19] O. Onireti, A. Mohamed, H. Pervaiz, and M. Imran, “Analytical approach to base station sleep mode power consumption and sleep depth,” in IEEE 28th Annu. Int. Symp. on Pers., Indoor, and Mobile Radio Commun. (PIMRC). Oct., 2017, pp. 1–7.
  * [20] J. P. Hespanha, “Modelling and analysis of stochastic hybrid systems,” IEE Proc.-Control Theory Appl., vol. 153, no. 5, pp. 520–535, Sep. 2006.


11institutetext: Purdue School of Industrial Engineering, West Lafayette, IN
47906, USA 11email<EMAIL_ADDRESS>22institutetext: Sardar Patel Institute
of Technology, Mumbai, 400058, India 22email<EMAIL_ADDRESS>33institutetext: D.J Sanghvi College of Engineering, Mumbai, 400056, India
33email<EMAIL_ADDRESS>

# Machine Learning to Estimate Gross Loss of Jewelry for Wax Patterns

Mihir Jain 11 Kashish Jain 22 Sandip Mane 33

###### Abstract

In mass manufacturing of jewellery, the gross loss is estimated before
manufacturing to calculate the wax weight of the pattern that would be
investment casted to make multiple identical pieces of jewellery. Machine
learning is a technology that is a part of AI which helps create a model with
decision-making capabilities based on a large set of user-defined data. In
this paper, the authors found a way to use Machine Learning in the jewellery
industry to estimate this crucial Gross Loss. Choosing a small data set of
manufactured rings and via regression analysis, it was found out that there is
a potential of reducing the error in estimation from ±2-3 to ±0.5- using ML
Algorithms from historic data and attributes collected from the CAD file
during the design phase itself. To evaluate the approach’s viability,
additional study must be undertaken with a larger data set.

###### Keywords:

CAD Gross Loss Jewelry Machine Learning Wax Model.

##  1 Introduction

Loss is an inevitable component of manufacturing. In the manufacturing of
jewellery from precious metals, accounting and calculating the losses is a
very crucial. Gross Loss of jewellery is the total metal loss during its
manufacturing. Loss in metal happens during casting, filing, polishing,
setting and at almost every stage. Even though most of this lost metal is
recovered and refined in the refinery to get a recovery of 92%, on average,
these losses are extremely crucial not to be accounted for.

The loss on each piece of jewellery varies, based on various factors.
Estimating this gross loss beforehand was very crucial for the manufacturing
of that jewellery. This estimated gross loss was used for while pulling wax
patterns during the process of injection moulding [1]. Jewelry made from the
heavier wax piece will have surplus metal that must be filed down and
recovered later, which is a waste of time and materials because only some of
the metal will be recovered. Therefore, estimating the total loss provides a
general estimate of the wax weight and can be used as a guide for how each
procedure should be carried out.

In a production process, a step wise loss of each of step of manufacturing is
collected. This is done by weighing the jewelry after each step. Hence after
the jewelry has been manufactured it can assess the final data of gross loss
that the company bore. Total recovery that was done was also considered, and
added to the database.

This gross loss found out was further collected out of which a wide set of
databases is manufactured by an in-house engineer. Calculations based on
current trends are made where a few other variables are also taken into
consideration.

Variables like, weight of the final product, metal type (White Gold, Yellow
Gold, Pink Gold, Silver, Platinum and Palladium), cartage of metal (8k, 9k,
10k, 12k, 14k, 18k, 20k etc.), the customer for whom the jewelry is being
manufactured, the setting of diamond (whether the piece is handset or wax set)
and of course the type of jewelry it is (whether it is a ring, a pendant, an
earring, a bracelet or a bangle.)

Currently, the estimation comes with a variance of ±4-5%. Hence there is a
scope here by which, using the powerful tools of Machine Learning [2, 3, 4, 5]
we can consider the variable constants to find out the gross loss in jewelry.
These variable constants can most often than not be fetched directly from the
CAD files which are made way before the actual manufacturing process even
begins.

The aim of the paper is to estimate the gross loss of jewelry at the CAD level
with greater and repeatable accuracy using machine learning algorithms. This
paper will systematically narrow down the variables responsible for gross loss
of jewelry during its manufacturing, create a machine learning model that
predicts the final gross loss based on the data collected from the CAD file
generated before manufacturing and ensure greater accuracy of the model as
compared to the traditional methods of estimating loss.

##  2 Methodology

As the project is a proof of concept, it only takes into account 26 rings as a
sample size. This project will only use information from the last several
months of production for all ring kinds for which CAD files were available
(developed in Rhino 3D [6]) and for which the company knew the associated
gross losses. It is important to highlight that only information that could be
shown publicly has been included in this report. There were notably three
stages to the project’s execution.

###  2.1 Creating the Dataset

The first phase comprised of selecting all possible attribute of the rings
from the CAD file and listing them down with their corresponding values on an
excel sheet. This data was paired with its corresponding historic gross loss.

Table 1: Parameters of the Dataset # | Attribute | Datatype
---|---|---
1 | Volume | mm3
2 | Surface Area | mm2
3 | Metal | Karat-metal
4 | Weight/ Piece (Estimated) | gm
5 | Total Lot Quantity | integer
6 | Total Weight of Lot | gm
7 | Inner Diameter | mm
8 | Outer Diameter | mm
9 | Minimum Shank Thickness | mm
10 | Maximum Shank Thickness | mm
11 | Minimum Shank Width | mm
12 | Maximum Shank Width | mm
13 | Total Height | mm
14 | Top Height | mm
15 | Number of Components | integer
16 | Number of Rings | integer
17 | Tone | 1/2/3
18 | True Miracle | binary
19 | No. of True Miracle | integer
20 | Diamond – Handset and Wax Set | integer
21 | Filigree | binary
22 | J Back | binary
23 | Gallery | binary
24 | Fake Beads | integer
25 | Plating | binary

###  2.2 Preparation of Data

The compiled data was obtained from the CAD files. This data had irrelevant
parameters that are currently unknown but will be filtered out through
processing. The reason why all possible data was collected was to avoid any
human generated discrepancies in the very first stage of the project. Even
though 26 is a small number for a machine learning algorithm, its
corresponding volume would still suffice to give us the proof of concept
required to carry on with the project. But in an ideal situation, the number
of rows should be 4x the number of columns. So, it can safely consider that
the results obtained will only improve when the amount of data of rings
increases. Before testing algorithms, the data was checked for any missing
values and such values were designated a weighted average value. Feature
scaling was done using Standard Scaling which standardize features by removing
the mean and scaling to unit variance [7]. The standard score of a sample x is
calculated as:

$z=(x-u)/s$ (1)

where u is the mean of the training samples, and s is the standard deviation
of the training samples.

###  2.3 Models

The above-mentioned data was split tested, and the mean value graph of each
feature was derived. The data set was split as 80-20, where 80 was the train
dataset and 20 was the test dataset. After the data set is split, feature
scaling was done to make the comparison easier. This data was the processed
through various algorithms.

  1. 1.

Linear Regression Algorithm

Figure 1: Architecture of Linear Regression Algorithm

    1. (a)

Linear Regression [8] is a supervised learning-based machine learning
technique. One of its functions is to carry out a regression analysis. Through
the use of independent variables, the regression model can predict a desired
outcome. Its primary function is to investigate causal links between factors
and predicting.

    2. (b)

Predicting a value (y) of a dependent variable (x) from known values (y) of
independent variables (x) is the job of linear regression (x). Accordingly,
this method of regression establishes a linear connection between the input
variable (x) and the output variable (y) (output). Linear Regression perfectly
describes this method [9]. Hypothesis function for Linear Regression :

$y=\theta_{1}+\theta_{2}.x$ (2)

    3. (c)

When training the model – it fits the best line to predict the value of y for
a given value of x. The model gets the best regression fit line by finding the
best $\theta$1 and $\theta$2 values where $\theta$1 is the intercept and
$\theta$2 is the coefficient of x.

    4. (d)

The best fit line is obtained by locating the optimal values of $\theta$1 and
$\theta$2. When our model is used for prediction, it will give us y as a
function of x.

  2. 2.

Random Forest Regression

Figure 2: Architecture of Random Forest Algorithm

    1. (a)

Random Forest Regression [10] algorithm is an example of a supervised learning
algorithm that use the ensemble learning approach of regression. By combining
the results of several different machine learning algorithms, an ensemble
learning method can produce a more precise forecast than any one of them could
on its own.

    2. (b)

Random Forest relies on the ”wisdom of the crowds” principle, which states
that a large number of independent models working in concert can achieve
better results than any of their parts working alone.

    3. (c)

This is owing to the fact that the trees buffer one another from their
particular errors. Since a random forest is completely random, there is no
communication between the trees that make up the forest. Random forest is an
estimator technique that takes the outputs of multiple decision trees,
compiles them, and then generates the ideal answer for the given situation.
[11]

  3. 3.

Decision Tree Regression

Figure 3: Architecture of ecision Tree Algorithm

    1. (a)

Decision tree [12] builds regression or classification models in the form of a
tree structure. It breaks down a dataset into smaller and smaller subsets
while at the same time an associated decision tree is incrementally developed.
The result is a tree with decision nodes and leaf nodes.

    2. (b)

There are three distinct sorts of nodes in this regression tree. The Root Node
is the primary node, representing the whole sample and potentially being
subdivided into further nodes. Features of a dataset are represented by
Interior Nodes, while decision rules are shown by Branches. In the end, the
result is shown by the Leaf Nodes. If you have an issue that requires a
choice, this algorithm is excellent. [13]

    3. (c)

A single data point is processed all the way down to the leaf node by asking
and answering True/False queries. Ultimately, the dependent variable value in
each leaf node is averaged to arrive at a final forecast. The Tree is able to
provide an accurate prediction after going through several rounds.

    4. (d)

The benefits of using decision trees include their simplicity, the fact that
they require less data cleaning, the fact that non-linearity has no effect on
the performance of the model, and the fact that the number of hyper-parameters
to be set is practically zero.

  4. 4.

K-Nearest Neighbors Regression

Figure 4: Example of KNN Algorithm

    1. (a)

K-Nearest Neighbors [14] is an easy-to-implement method that remembers past
examples and makes a numerical prediction based on a similarity measure (e.g.,
distance functions). KNN is a non-parametric technique that has been utilised
for statistical estimates and pattern recognition since the early 1970s.

    2. (b)

We need to use cross-validation to choose K. Unlike classification we cannot
accuracy use as metric, since our predictions will almost never exactly match
the true response variable values. Therefore in the context of KNN regression
we will use root mean square prediction error (RMSPE) instead. The
mathematical formula for calculating RMSPE is:

$RMSPE=\sqrt{\frac{\sum_{1}^{n}\left(y_{i}-\hat{y_{i}}\right)^{2}}{n}}$ (3)

Where n is the number of observations, yi is the observed value for the ith
observation, and ŷ\textsubscript{i} is the predicted value for the ith
observation.

    3. (c)

To put it another way, for each observation in our test (or validation) set,
we calculate the squared difference between the anticipated and true response
value, average over observations, and then square root. Since differences can
be either positive or negative—that is, we can over- or under-estimate the
genuine response value—we utilise the squared difference rather than merely
the difference. [15]

##  3 Results

As we’ve seen, we have put our datasets through their paces with a wide range
of ML algorithms. The Mean Absolute Error (MAE) represents the error that was
introduced into our findings. The formula for the same is

$MAE=\frac{\sum_{i=1}^{n}|y_{i}-x_{i}|)}{n}$ (4)

Where, yi is the prediction, xi is the true value and n specifies the total
number of data points.

Table 2: Mean Absolute Error Of Different Methods Method | Mean Absolute Error
---|---
Linear Regression | 0.56
Random Forest Regressor | 1.72
Decision Tree Regressor | 1.49
K-Nearest Neighbour Regressor | 2.02

It was observed that all 4 algorithms performed well considering how small the
data set was. All algorithms gave promising results with Linear Regression
lending the lowest MAE. Though with increasing data set, it would be wise to
consider all the remaining models as well. The scores will only improve as the
data set increases.

##  4 Conclusion

The results show us that the gross loss can be predicted to an error margin of
±0.5. The proof of concept needed to be derived from the results was
sufficient to take to the company to act on it. Each of the 4 algorithms have
potential, Linear Regression being the most promising one so far. Further
testing needs to be done by increasing the number of data set and even
expanding to different category of jewelry. The implementation of this
innovation in the field of jewelry manufacturing would be a big undertaking
and would be a time consuming and labour-intensive processes but one which
would bare fruitful results.

## References

  * [1] F.R. Sias.  Lost-wax Casting: Old, New, and Inexpensive Methods.  Woodsmere Press, 2005.
  * [2] S.J. Raudys and A.K. Jain.  Small sample size effects in statistical pattern recognition: recommendations for practitioners.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 13(3):252–264, 1991.
  * [3] Christopher Bishop.  Pattern Recognition and Machine Learning, volume 16, pages 140–155.  01 2006.
  * [4] Pedro Domingos.  A few useful things to know about machine learning.  Commun. ACM, 55(10):78–87, oct 2012.
  * [5] Ian H. Witten, Eibe Frank, and Mark A. Hall.  Data Mining: Practical Machine Learning Tools and Techniques.  Morgan Kaufmann Publishers Inc., San Francisco, CA, USA, 3rd edition, 2011\.
  * [6] Robert McNeel et al.  Rhinoceros 3d, version 6.0.  Robert McNeel & Associates, Seattle, WA, 2010.
  * [7] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay.  Scikit-learn: Machine learning in Python.  Journal of Machine Learning Research, 12:2825–2830, 2011.
  * [8] David Freedman.  Statistical Models: Theory and Practice.  01 2005.
  * [9] Khushbu Kumari and Suniti Yadav.  Linear regression analysis study.  Journal of the Practice of Cardiovascular Sciences, 4:33, 01 2018\.
  * [10] Tin Kam Ho.  Random decision forests.  In Proceedings of the Third International Conference on Document Analysis and Recognition (Volume 1) - Volume 1, ICDAR ’95, page 278, USA, 1995\. IEEE Computer Society.
  * [11] Ashwin Raj.  A quick and dirty guide to random forest regression, Jun 2021.
  * [12] Shai Shalev-Shwartz and Shai Ben-David.  Decision Trees, page 212–218.  Cambridge University Press, 2014.
  * [13] Gurucharan M K.  Machine learning basics: Decision tree regression, Jul 2020.
  * [14] Ruidi Chen and Ioannis Ch. Paschalidis.  Selecting Optimal Decisions via Distributionally Robust Nearest-Neighbor Regression.  Curran Associates Inc., Red Hook, NY, USA, 2019.
  * [15] Kevin Beyer, Jonathan Goldstein, Raghu Ramakrishnan, and Uri Shaft.  When is “nearest neighbor” meaningful?  In Catriel Beeri and Peter Buneman, editors, Database Theory — ICDT’99, pages 217–235, Berlin, Heidelberg, 1999. Springer Berlin Heidelberg.


11institutetext:  aKavli IPMU (WPI), UTIAS, The University of Tokyo, Kashiwa,
Chiba 277-8583, Japan, bDepartment of Physics, University of Tokyo, 7-3-1,
Hongo, Bunkyo-ku, Tokyo 113-0033, Japan, cOkayama University, 3-1-1,
Tsushimanaka, Kita-ku, Okayama City, Okayama, 700-8530, Japan, dSaitama
University, 255 Shimookubo, Sakura-ku, Saitama, 338-8570, Japan.
$\dagger$Corresponding author<EMAIL_ADDRESS>
# Modelling TES non-linearity induced by a rotating HWP in a CMB polarimeter

T. Ghignaa, $\dagger$ T. Matsumuraa Y. Sakuraia R. Takakub K. Komatsuc S.
Sugiyamad Y. Hoshinod N. Katayamaa

###### Abstract

Most upcoming CMB experiments are planning to deploy between a few thousand
and a few hundred thousand TES bolometers in order to drastically increase
sensitivity and unveil the B-mode signal. Differential systematic effects and
$1/f$ noise are two of the challenges that need to be overcome in order to
achieve this result. In recent years, rotating Half-Wave Plates have become
increasingly more popular as a solution to mitigate these effects, especially
for those experiments that are targeting the largest angular scales. However,
other effects may appear when a rotating HWP is being employed. In this paper
we focus on HWP synchronous signals, which are due to intensity to
polarization leakage induced by a rotating cryogenic multi-layer sapphire HWP
employed as the first optical element of the telescope system. We use LiteBIRD
LFT as a case study and we analyze the interaction between these spurious
signals and TES bolometers, to determine whether this signal can contaminate
the bolometer response. We present the results of simulations for a few
different TES model assumptions and different spurious signal amplitudes.
Modelling these effects is fundamental to find what leakage level can be
tolerated and minimize non-linearity effects of the bolometer response.

###### Keywords:

CMB, Bolometers, Transition Edge Sensors, Half-Wave Plates

††journal: Journal of Low Temperature Physics

##  1 Introduction

Probing the primordial CMB B-mode signal is one of the main objectives of
modern cosmology. While still an open question, a very rapid inflationary
expansion of the Universe is expected from current theories. The best
available tool to probe inflation is a measurement of the polarized CMB signal
zaldarriaga_seljak97 ; kamionkowski97 . On large angular scales ($\ell\lesssim
200$) the polarized B-mode signal is expected to be dominated by a primordial
inflationary component (if the tensor-to-scalar ratio $r\gtrsim 0.01$, or even
$r\gtrsim 0.001$ for $\ell<10$)111Galactic foregrounds and gravitational
lensing are two major challenges for a B-mode measurement. While a detailed
discussion is beyond the scope of this paper, we briefly mention them here
because their presence can make non-linear effects even more relevant.. Two of
the main challenges for CMB experiments are atmospheric fluctuations (for
ground-based experiments) and receiver stability on long time-scales ($1/f$
noise and gain fluctuations). Therefore, several present and future CMB
experiments (CLASS Dahal2021 , Simons Array Suzuki2016 , Simons Observatory
simonsObs2019 , LiteBIRD Hazumi2020 ) are employing (or planning to)
polarization modulators, such as continuously-rotating Half-Wave Plates (HWP).
A polarization modulator up-converts the polarized sky signal to a frequency
range where the noise spectrum is expected to be purely white Kusaka2014 . In
a nutshell, if the HWP rotates at a frequency $f_{hwp}$, the polarized signal
is modulated at $f_{mod}=4f_{hwp}$, therefore by properly choosing $f_{hwp}$
it is possible to shift the signal of interest above the $1/f$ knee-frequency.

LiteBIRD is one of the experiments that will adopt this strategy. The LiteBIRD
space mission consists of 3 telescopes each equipped with a HWP, as the first
optical element, and a cryogenic focal plane ($\sim 100$ mK) populated with
polarization sensitive Transition-Edge Sensor (TES) bolometers. At Kavli IPMU
we are in charge of developing the HWP and the rotation mechanism for LFT. The
current baseline consists of a multi-layer sapphire HWP with anti-reflection
layers laser-machined directly on the HWP surfaces Komatsu2019 ; Takaku2021 .
Even though HWPs are effective in rejecting $1/f$ noise, as well as removing
some systematic effects arising from intrinsic differences among detectors
(e.g. bandpass and beam mismatch, gain mismatch, etc.), HWP non-idealities can
introduce other systematic effects. For example, in Takakura2017 ; Didier2017
(POLARBEAR and EBEX collaborations) the authors identified large222Amplitudes
of the order of hundreds of mK. For comparison the CMB solar dipole is
$\sim\pm 3.362$ mK. signals due to instrumental intensity-to-polarization
leakage. In both cases the signal is attributed primarily to optical elements
before the HWP, and it is modulated by the HWP at $4f_{hwp}$.333Any optical
element located before the HWP generates instrumental polarization that is
going to be modulated by the HWP at the same frequency of the sky signal.
Having the HWP as the first optical element allows for a clear separation of
instrumental and celestial polarization. In both analysis the authors argue
that this signal affects the TES response by driving it into a non-linear
regime.

In the case of LiteBIRD and Simons Observatory the HWP is going to be the
first optical element (if we neglect windows or thermal filters), therefore we
should not expect the same $4f_{hwp}$ leakage level, however HWP non-
idealities or non-orthogonal incidence can create similar effects.

In this paper we address the problem of TES non-linearity and highlight how
this can couple to HWP synchronous signals (HWPSS). We use LiteBIRD LFT as a
case study. First, we define TES non-linearity and HWP intensity-to-
polarization leakage.

(a)

(b)

Figure 1: Example of the input maps in pW at 100 GHz. The left plot shows the
total intensity map including CMB anisotropies, foreground (generated with
pysm) the CMB monopole and dipole and the telescope loading for a LiteBIRD LFT
100 GHz detector according to Westbrook2021 ; Hasebe2021 . The right plot
shows the polarization map $\sqrt{Q^{2}+U^{2}}$. The maps are generated in CMB
units and then converted to input power. (Color figure online.)

(a)

(b)

Figure 2: Left: 1 hour of noiseless time-ordered data without IP leakage
(ideal case) generated from the map in Figure 1 according to LiteBIRD scan
strategy litebird_ptep  for a detector at the telescope boresight. Right:
time-ordered data after processing the data on the left with the TES detector
simulator. The data contains both total intensity and polarization. The
presence of the dipole and galactic plane can be easily seen in the time-
stream. (Color figure online.)

We present a code developed in python to realistically simulate the detector
response to the sky signal. Finally, we present preliminary results showing
the impact of spurious signals (such as HWPSS) in the detector response. We
compare these results to non-linearity estimates from POLARBEAR data
Takakura2017 . We believe that the tool presented can be useful to forecast
future experiment capabilities, define instrument requirements, inform the
instrument design and develop methodologies to analyze the data and mitigate
systematic effects. This is vital for the level of accuracy required by future
experiments in particular for a space mission like LiteBIRD, which is aiming
to achieve an accuracy in terms of tensor-to-scalar ratio of
$\sigma_{r}\lesssim 0.001$. A detailed study of LiteBIRD requirements can be
found in litebird_ptep .

##  2 TES and HWP model

### Transition-Edge Sensor:

According to Irwin2005  we can define the behaviour of a DC voltage-biased
bolometer by solving the system of coupled differential equations:

$\displaystyle L\frac{dI}{dt}$ $\displaystyle=$ $\displaystyle
V_{b}-IR_{tes}-IR_{sh}$ (1) $\displaystyle C\frac{dT}{dt}$ $\displaystyle=$
$\displaystyle-P_{b}+I^{2}R_{tes}+P_{opt}.$ (2)

It is commonly understood that, in the small signal ($\delta
P_{opt}\ll\bar{P}_{opt}$) and high loop-gain limits $\mathcal{L}\ll 1$, the
detector current responsivity reduces to $S_{I}\sim-1/V_{b}$ (for a ”slow”
signal: $\omega\rightarrow 0$).

In Takakura2017 ; TakakuraPHD  the authors present a non-linearity model where
the detector responsivity is loosely dependent from the input power variations
$\delta P_{opt}$. This takes the following form for the detector current:

$\delta I(t)=[S_{I}+S^{\prime}_{I}\delta P_{opt}(t)]\delta
P_{opt}(t-\tau^{\prime}\delta P_{opt}(t)),$ (3)

where $S^{\prime}_{I}$ and $\tau^{\prime}$ are non-linear terms. In the
following, after briefly describing the HWP, we will discuss how the non-
linear term couples to the HWP creating spurious signals.

### Half-Wave Plate:

The effect of a rotating HWP on the incoming sky signal is defined in terms of
its Mueller matrix $\Gamma_{hwp}$ Komatsu2019 . For an ideal HWP
$\Gamma_{hwp}$ is a diagonal matrix of the form $diag(1,1,-1,-1)$. However,
for a real HWP the off-diagonal elements of the matrix $\Gamma_{hwp}$ can be
non-zero, and cause mixing between the $I$, $Q$ and $U$ components of the
input signal. In particular, in the context of this paper,

(a)

(b)

(c)

(d)

Figure 3: Top: 1% IP leakage. Left: NSD (noise spectral density) of the input
TOD in units of power (1 hour). Right: NSD of the TOD after processing with
the TES simulator. Bottom: Same as top for 10% IP leakage. These simulations
are noiseless, therefore the spectrum is dominated by the foreground signal.
(Color figure online.)

we are interested in the $\gamma_{QI}$ and $\gamma_{UI}$ components that can
cause intensity-to-polarization leakage (IP). From Komatsu2019  we find that
non-zero values can result in a $2f_{hwp}$ signal of the form:
$\frac{1}{2}\sqrt{\gamma^{2}_{QI}+\gamma^{2}_{UI}}I_{in}\cos(2\omega_{hwp}t-2\arctan{\frac{\gamma_{UI}}{\gamma_{QI}}})$
Takaku2021 . Other effects can also appear, if we take into consideration the
position of the HWP in the optical system, but for the scope of this paper we
limit the discussion to the HWP intensity-to-polarization leakage. However,
this analysis can be easily generalized to other cases.

##  3 Simulation

As discussed in Section 2, if the TES responsivity depends loosely on the
incoming signal $\delta P_{opt}$ the bolometer response becomes non-linear.
From Equation 3 we can easily see that a signal at a generic frequency
$\omega$ can partially leak to $2\omega$.444From simple trigonometry, if
$\delta P_{opt}\sim\sin{\omega t}$, then the term $(\delta
P_{opt})^{2}\sim\frac{1-\cos{2\omega t}}{2}$. Hence, the $2f_{hwp}$ IP signal
discussed in Section 2 can partially leak to $4f_{hwp}$, where the polarized
signal is expected. In order to perform these simulations and address the
impact of non-linearity on cosmological data we have developed a $python$
module555https://github.com/tomma90/tessimdc that can be plugged into a more
standard time-order data (TOD) simulator. First, the input sky maps are
converted to units of power by integrating the sky emission ($Jy/sr$) over the
beam function $B(\Omega,\nu)$ and the band-pass function $G(\nu)$. An example
is shown in Figure 1 for LiteBIRD LFT 100 GHz. We use the
$pysm$666https://pysm3.readthedocs.io/en/latest/ package to generate the sky
model. We then add the CMB dipole and the expected optical loading (for a
LiteBIRD LFT 100 GHz detector) due to the CMB monopole and telescope loading:
$\bar{P}_{opt}=0.3061$ pW Westbrook2021 ; Hasebe2021 .

The second step consists in generating the telescope pointing according to
LiteBIRD scan strategy Hazumi2020  which is used to generate the TOD in units
of power. At this point in the simulation, we can add a HWP synchronous
signal.

Figure 4: TES detector non-linearity. We simulated representative sinusoidal
signals of frequencies $f_{sig}=$ 0.1, 2 and 4 Hz and amplitudes between 0.1%
and 10% of the expected average optical loading $P_{opt}$ to test non-
linearity according to the formula in Equation 3. We tested 2 cases: a TES
with log-responsivity $\alpha_{0.5\Omega}$=100 and $\alpha_{0.5\Omega}$=1000
at $R=0.5\times R_{n}$. For each case we biased the TES at $R=0.5\times R_{n}$
and $R=0.75\times R_{n}$ and studied its response. It is clear that a narrower
transition and a higher loop-gain lead to a more linear device. (Color figure
online.)

At present we can only fix the amplitude of the signal. We are planning to
implement the possibility to input a measured or simulated HWP Mueller matrix
to allow for more realism to the simulation. The TOD is then fed to the TES
simulator as the input optical power $P_{opt}$ in Equation 2. An example for a
1-hour scan is shown in the left plot of Figure 2.

The TES response simulator consists of a module based on the Runge-Kutta 4th
order method to solve Equations 1 and 2. We have made some assumptions to set
the parameters according to LiteBIRD specifications. The thermal conductance
$G$ and heat capacity $C$ are defined assuming a normal time-constant
$\tau_{0}=C/G\sim 30$ ms and a saturation power
$P_{sat}=2.5\times\bar{P}_{opt}$ (average optical power defined by the CMB
monopole and telescope loading).

Another key assumption is the TES resistance. In the simulations carried out
for this paper, we have assumed a simple $arctan$ model where the resistance
depends only on the TES temperature. We defined the normal resistance
$R_{n}=1\Omega$ and the log-responsivity ($\alpha=\frac{dlogR}{dlogT}$)
$\alpha_{0.5\Omega}=100$ at $0.5\times R_{n}$ (conservative value). In the
future we will explore other models that do not neglect the dependence of the
resistance from the TES current (two-fluids, RSJ, etc.) One shortcoming of
this procedure is the need for a high sampling rate compared to more
traditional TOD simulator. This is due to the need to solve iteratively the
differential equations at a rate faster than the detector effective time-
constant ($\tau_{eff}\sim\tau_{0}/(\mathcal{L}+1)$) in order to avoid
numerical errors. For high-loop gain conditions ($\mathcal{L}\sim 10$), this
translates to a sampling rate $\gtrsim$ 1 kHz.

In the right plot of Figure 2 we can see the result of the TES simulator for
an ideal case without HWPSS. In Figure 3 we show the impact of a $2f_{hwp}$
($\omega_{hwp}=46$ rpm) HWPSS component of $1\%$ and $10\%$ amplitudes
(relative to the expected optical loading $\bar{P}_{opt}$)777The data used in
this work are noiseless to focus on the detector response characterization.. A
larger $2f_{hwp}$ signal clearly results into an enhancement of the $4f_{hwp}$
component in the NSD (noise spectral density), which is due to the non-linear
term of Equation 3 that grows with the amplitude of the signal
($S^{\prime}_{I}\delta P^{2}$).

Finally, we test the non-linearity model in Equation 3 using representative
signals (for this simulation we simplified the input compared to Figure 3, by
including only the monopole - $P_{opt}$ \- and a signal of frequency $f_{sig}$
and relative amplitude $\delta P/P_{opt}$) of frequency $f_{sig}$ between 0.1
Hz and 4 Hz and amplitudes between 0.1% and 10% of the expected optical
loading for two transition cases: a log-responsivity $\alpha_{0.5\Omega}$=100
and $\alpha_{0.5\Omega}$=1000. For all cases we study the response when the
detector is biased at $R=0.5\times R_{n}$ and $R=0.75\times R_{n}$. After
processing the input signal with the TES simulator we convert the TES current
$\delta I$ to input power $\delta P^{\prime}=\delta I/S_{I}$, in the high
loop-gain limit. Afterwards we compare the input signal $\delta P$ and the
calibrated output $\delta P^{\prime}$ to determine the non-linear term
$S^{\prime}_{I}\delta P$ for each case: $S^{\prime}_{I}\delta P/S_{I}=(\delta
P^{\prime}-\delta P)/\delta P$. In Figure 4 we show a summary of the results.
By isolating the term $S^{\prime}_{I}\delta P$ we find that this appears to
depend on the TES model and bias conditions. In fact, the cases with steeper
and narrower transition ($\alpha_{0.5\Omega}$=1000) and lower operating
resistance (higher loop-gain) show reduced non-linearity levels.

##  4 Conclusions

We believe that the tool developed and presented in this paper is going to be
very valuable going forward with the definition of the mission, both in terms
of informing the instrument development and the data analysis. This is
particularly true for the requirements of a space mission. We have given a
preliminary assessment of the TES non-linearity model developed and tested on
POLARBEAR data in TakakuraPHD ; Takakura2017 . We found indications of the
input power-dependent non-linear component. While the responsivity term seems
to mimic well the model, we need a more in depth analysis of the non-linearity
model to address the time-constant term. A more thorough analysis with
propagation to map-making and cosmological analysis is needed to understand
the full impact of this effect on the scientific output (e.g. tensor-to-scalar
ratio). The results888The data that support the findings of this study are
available from the corresponding author upon reasonable request. presented in
this paper are a starting point for a more in-depth analysis. More data and
laboratory tests are needed to define realistic detector and instrument models
for LiteBIRD.

However, we can notice that in Takakura2017  the level of non-linearity of
POLARBEAR data has been estimated to be in the range 0.3-0.8% (varying among
detector wafers). This is in agreement with the values found in our analysis
for the cases with $\alpha_{0.5\Omega}=1000$ in Figure 4, where we have
determined $S^{\prime}_{I}\delta P/S_{I}$ to be in the range 0.2-0.9%
(depending on the bias conditions). In our analysis we have also studied the
case for $\alpha_{0.5\Omega}=100$ which is regarded as a pessimistic scenario,
and in fact leads to larger non-linearity levels (as high as $\sim 9$%
according to Figure 4). This result clearly shows that a narrower
superconductive transition improves the linearity of the device.

In conclusion, the absence of a turbulent atmosphere that affects ground-based
experiments will allow LiteBIRD (or any other space mission) to achieve an
unprecedented sensitivity. We can foresee that more stable observation
conditions will reduce the impact of $1/f$ noise and gain fluctuations due to
long time-scale variations of the optical loading. Hence, more subtle non-
linear effects could become more dominant, which will require more realistic
models of the detectors and read-out systems to forecast and mitigate them.

###### Acknowledgements.

We thank all LiteBIRD collaborators for support and help. In particular Juan
Macías-Pérez and Satoru Takakura for usuful comments and feedback on the
manuscript. This work was supported by JSPS KAKENHI Grant Numbers 22K14054,
18KK0083 and 19K14732. Kavli IPMU is supported by World Premier International
Research Center Initiative (WPI), MEXT, Japan.

## References

  * (1) M. Zaldarriaga and U. Seljak. Phys. Rev. D 55, 4, 1830-1840 (1997). DOI: 10.1103/PhysRevD.55.1830.
  * (2) M. Kamionkowski, A. Kosowsky, and A. Stebbins. Phys. Rev. Lett. 78, 11, 2058-2061 (1997). DOI: 10.1103/PhysRevLett.78.2058.
  * (3) Dahal S., Appel J. W., Datta R., Brewer M. K., Ali A., Bennett C. L., Bustos R., et al., 2021, arXiv, arXiv:2107.08022.
  * (4) Suzuki A., Ade P., Akiba Y., Aleman C., Arnold K., Baccigalupi C., Barch B., et al., 2016, JLTP, 184, 805. doi:10.1007/s10909-015-1425-4.
  * (5) Ade P., Aguirre J., Ahmed Z., Aiola S., Ali A., Alonso D., Alvarez M. A., et al., 2019, JCAP, 2019, 056. doi:10.1088/1475-7516/2019/02/056.
  * (6) Hazumi M., Ade P. A. R., Adler A., Allys E., Arnold K., Auguste D., Aumont J., et al., 2020, SPIE, 11443, 114432F. doi:10.1117/12.2563050.
  * (7) Kusaka A., Essinger-Hileman T., Appel J. W., Gallardo P., Irwin K. D., Jarosik N., Nolta M. R., et al., 2014, RScI, 85, 039901. doi:10.1063/1.4867655.
  * (8) Komatsu K., Matsumura T., Imada H., Ishino H., Katayama N., Sakurai Y., 2019, JATIS, 5, 044008. doi:10.1117/1.JATIS.5.4.044008.
  * (9) R. Takaku et al, , J. Low Temp. Phys., This special issue (2021).
  * (10) Takakura S., Aguilar M., Akiba Y., Arnold K., Baccigalupi C., Barron D., Beckman S., et al., JCAP, 05, 2017, 008. doi:10.1088/1475-7516/2017/05/008.
  * (11) Didier J., Miller A. D., Araujo D., Aubin F., Geach C., Johnson B., Korotkov A., et al., 2019, ApJ, 876, 54. doi:10.3847/1538-4357/ab0f36.
  * (12) T. Hasebe et al, , J. Low Temp. Phys., This special issue (2021).
  * (13) Westbrook B., Raum C., Beckman S., Lee A. T., Farias N., Sasse T., Suzuki A., et al., 2021, arXiv, arXiv:2101.05306.
  * (14) LiteBIRD Collaboration: Allys, E., Arnold, K., Aumont, J., Aurlien, R., Azzoni, S. et al., 2022, arXiv, arXiv:2202.02773.
  * (15) K. Irwin, G. Hilton, Ch: Transition-Edge Sensors. In: Cryogenic Particle Detection. Topics in Applied Physics, DOI: 10.1007/10933596_3.
  * (16) S. 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, DOI: 10.18910/67054.
  * (17) J. C. Butcher, Numerical Methods for Ordinary Differential Equations, 2008, DOI: 10.1002/9780470753767.