# Average AoI Minimization for Energy Harvesting Relay-aided Status Update
Network Using Deep Reinforcement Learning

Sin-Yu Huang and Kuang-Hao (Stanley) Liu S-Y. Huang is with the Institute of
Computer and Communication Engineering, National Cheng Kung University,
Tainan, Taiwan 701. E-mail<EMAIL_ADDRESS>Liu is with the
Institute of Communications Engineering, National Tsing Hua University,
Hsinchu, Taiwan 300044. E-mail<EMAIL_ADDRESS>
###### Abstract

A dual-hop status update system aided by energy-harvesting (EH) relays with
finite data and energy buffers is studied in this work. To achieve timely
status updates, the best relays should be selected to minimize the average age
of information (AoI), which is a recently proposed metric to evaluate
information freshness. The average AoI minimization can be formulated as a
Markov decision process (MDP), but the state space for capturing channel and
buffer evolution grows exponentially with the number of relays, leading to
high solution complexity. We propose a relay selection (RS) scheme based on
deep reinforcement learning (DRL) according to the instantaneous channel
packet freshness and buffer information of each relay. Simulation results show
a significant improvement of the proposed DRL-based RS scheme over state-of-
art approaches.

###### Index Terms:

Age of information, buffer-aided relaying, energy harvesting, relay selection,
status update.

##  I Introduction

Emerging Internet of Things (IoT) stimulates the need for timely status
updates. For example, health and meteorological data need to be updated using
the data collected by remote sensors in a timely manner [1, 2]. To this end,
the communication network should be designed to ensure timely status update,
which can not be directly achieved through optimizing traditional metrics,
such as outage probability, throughput, and queueing delay. Recently, a new
metric, namely the _Age of Information_ (AoI) [3], is proposed to evaluate the
freshness of the received data. Existing studies on AoI focus on wireline
networks with fruitful theoretical insights established based on queueing
theories [4] and recent attention have been turned to wireless networks [5].
Clearly, timely delivery of status packets is challenging in wireless networks
due to time-varying wireless channel conditions.

In many IoT applications, the information sink node might be distant from the
sensing devices or their direct links may be blocked by large obstacles that
make timely status updates difficult. To remedy, intermediate relays can be
deployed to extend the effective communication range of the sensing devices
with limited transmission power. Information relaying using cooperative relays
has been shown to achieve spatial diversity without the need for a large
antenna array [6], which is not practical for sensing devices with small form
factors. Some existing work has investigated how to enable status updates
using relays. In [7], the authors analyze the AoI of the wireless relay
network based on the theoretical tool of Markovian jump linear systems. A
static link scheduling policy, which assigns the transmission opportunity to
either the source or the relay nodes with a fixed probability, is considered
in their work. Since the transmission probability for AoI minimization does
not have a closed form, it is optimized using numerical search in [7]. In [8,
9], conventional relay selection (RS) schemes are applied to the dual-hop
status update network with a focus on the tradeoff between the average AoI and
the average energy cost. In the aforementioned work, relays are assumed to
have a permanent power supply, which requires frequent battery change or
costly power cabling. For IoT use cases, an appealing implementation choice is
to deploy EH relays that replenish energy from the source signal. [10]
explores the use of EH relay for real-time status updates in IoT by leveraging
the notion of simultaneous wireless information and power transfer (SWIPT).
However, their work is limited to a single-relay case and thus no spatial
diversity is exploited for status updates.

In light of the explosive demand for real-time sensing in IoT, this work
studies the AoI minimization problem in the dual-hop status update network
that has not been well addressed. In this context, efficient status updates
can benefit from selection diversity, namely choosing the best relay to
receive and transmit the status update packet, respectively. Unlike [8, 9]
where the considered RS schemes are based on channel status information (CSI)
without a specific design for ensuring information freshness, the RS problem
for AoI minimization is solved in this work by a model-free approach based on
deep reinforcement learning (DRL). The proposed DRL-based RS scheme makes the
RS decision by jointly considering CSI packet freshness and buffer evolution
relevant to the AoI. The advantage of the proposed DRL-based RS scheme is that
it bypasses the difficulty of handling a large state space encountered in the
Markov decision process (MDP) based formulation in order to arrive the optimal
RS decision.

The remainder of this paper is organized as follows. Sec. II illustrates the
system model and problem formulation. The proposed DRL-based RS scheme is
described in Sec. III. Sec. IV presents the numerical results. Sec. V
summarizes the key results. Reproducible research: All the simulation results
can be reproduced using Python code available at:
https://github.com/Sinyu104/Relay-Selection-based-on-DQN.

##  II System Model and Problem Formulation

We consider a dual-hop status update network where a source node $S$ sends
status update packets to a destination node $D$ with the aid of $K$
intermediate relay nodes $R_{1},\cdots,R_{K}$ using decode-and-forward (DF). A
status packet is generated in a slot with probability $\lambda$. The slot
length is chosen to allow at most one packet transmission. All nodes are
single-antenna half-duplex devices. The direct communication between $S$ and
$D$ is not available due to severe blockage. Therefore, $D$ can only receive
data via relays that are self-sustainable by replenishing energy from the
radio frequency (RF) signal sent from $S$. To this end, each relay has an
energy-harvesting (EH) circuit that converts the RF signal into DC with a
conversion efficiency $\eta\in(0,1)$ [11]. The received packet and energy is
stored in a data buffer and an energy buffer, respectively, both with a finite
storage [12]. The data buffer can accommodate at most one packet, because
backlogged packets in the data buffer is not helpful to status update [13].
Each packet stored in $R_{k}$ consists of a time stamp $\alpha_{k}$ recording
its generation time. The energy buffer has a size of $E_{\max}$ in the unit of
energy intervals. Here, an energy interval is equal to the energy consumption
for transmitting one packet by the relay. The instantaneous lengths of the
data and the energy buffers of relay $R_{k}$ are denoted as
$\delta_{k}\in\\{0,1\\}$ and $\epsilon_{k}\in\\{0,1,\cdots,E_{\max}\\}$,
respectively. The constant power transmitted by the source and the relays is
denoted as $P_{S}$ and $P_{R}$, respectively.

The transmission over each link experiences path-loss and multi-path fading.
Accordingly, the channel coefficients between $S-R_{k}$ and $R_{k}-D$ links,
denoted as $h_{k}$ and $g_{k}$ for $k=1,\cdots,K$, are modeled as independent
and identically distributed (i.i.d) zero-mean complex Gaussian random
variables with variances $d^{-2}_{h_{k}}$ and $d^{-2}_{g_{k}}$ with
$d_{h_{k}}$ and $d_{g_{k}}$ being the distance between $S$ and $R_{k}$ and
that between $R_{k}$ and $D$, respectively. The instantaneous signal-to-noise
ratio (SNR) of the $S-R_{k}$ link is symbolized as
$\gamma_{h_{k}}\triangleq|h_{k}|^{2}\frac{P_{S}}{N_{0}}$ where $N_{0}$ denotes
the noise power. Similarly, the instantaneous SNR of the $R_{k}-D$ link is
represented as $\gamma_{g_{k}}\triangleq|g_{k}|^{2}\frac{P_{R}}{N_{0}}$.
Denote $t$ the slot length, corresponding to the time duration for
transmitting one packet. The amount of energy harvested by $R_{k}$ in one slot
is given by

$\displaystyle E_{k}=\eta P_{S}\left|h_{k}\right|^{2}\cdot t.$ (1)

###  II-A Transmission Protocol

To leverage relays for status update, the following protocol is used. In each
time slot, the DRL agent at $S$ decides whether to activate an $S-R$ link or
an $R-D$ link. If an $S-R$ link is activated, $S$ transmits a packet to a
selected relay and the remaining relays replenish energy from the source
signal. The received packet by the relay will be stored in the data buffer
before transmission. If an $R-D$ link is activated, a selected relay will
forward a packet to $D$ while other relays remain idle.

In measuring information freshness at $D$, we use AoI as the metric [14]. For
slot $n$, the AoI denoted as $A(n)$ is the time elapsed since the generation
time of the last successfully received packet by $D$, i.e,
$A(n)=n-\hat{\alpha}(n)$, where $\hat{\alpha}(n)$ is the generation moment of
the latest received update packet by $D$ at slot $n$. The packet forwarded
from a relay might be outdated, i.e.,
$\hat{\alpha}(n^{\prime})<\hat{\alpha}(n)$ for $n^{\prime}<n$ in which case
AoI will be accumulated. Otherwise, AoI is reset to the time that elapsed
since the generation of the latest packet. In summary, the AoI at slot $n$ is
given by

$\displaystyle
A(n)=\begin{cases}n-\hat{\alpha}(n),&~{}\hat{\alpha}(n)>\hat{\alpha}(n^{\prime}),n>n^{\prime}\\\
\min\\{A(n-1)+1,A_{\max}\\},&~{}\text{otherwise}.\end{cases}$ (2)

where $A_{\max}$ is a prescribed upper bound. A packet received at $D$ with
the AoI exceeding $A_{\max}$ will be discarded. Also, the packet in the relay
buffer might be outdated when it is generated later than that of the latest
received packet by $D$. To measure the freshness of the buffered packet at
relay $R_{k}$, define the relative age as
$\alpha^{\prime}_{k}(n)\triangleq\alpha_{k}(n)-\hat{\alpha}(n)$ [13].

###  II-B Problem Formulation

The agent at $S$ aims to minimize the average AoI by selecting either a
reception relay denoted as $R_{r}(n)$ or a transmitting relay denoted as
$R_{t}(n)$ in slot $n$. For a time window of $N$ slots, the average AoI
minimization problem can be formulated as follows.

$\displaystyle(\text{P})$
$\displaystyle\quad\mathop{\min}\limits_{\boldsymbol{\mathcal{R}}}~{}\frac{1}{N}\sum_{n=1}^{N}A(n)$
s.t. $\displaystyle\quad\delta_{k}(n)\in\\{0,1\\},$ (3a) $\displaystyle\quad
0\leq\epsilon_{k}(n)\leq E_{\max},$ (3b)
$\displaystyle\quad\delta_{R_{t}}(n)=1~{}\text{and}~{}\epsilon_{R_{t}}(n)\geq
1,$ (3c)
$\displaystyle\quad\gamma_{R_{r}}(n)\geq\gamma_{\text{th}}~{}\text{and}~{}\gamma_{R_{t}}(n)\geq\gamma_{\text{th}},$
(3d) $\displaystyle\quad\alpha^{\prime}_{R_{t}}(n)>0.$ (3e)

for $k=1,\cdots,K$ and $n=1,\cdots,N$. In (P), $\boldsymbol{\mathcal{R}}$ is
an $N\times 1$ vector with each entry indicating the index of the selected
relay. Constraints (3a) and (3b) follow due to the size limit of the data and
energy buffers, respectively. Constraint (3c) ensures that the transmitting
relay has a non-empty data buffer and its available energy is sufficient for
transmitting a packet. The activation of an $S-R$ link and an $R-D$ link takes
place when their instantaneous SNRs exceeds a threshold $\gamma_{\text{th}}$,
as stated in (3d). Finally, constraints (3e) guarantees that an outdated
packet would not be sent to $D$. Problem (P) can be treated as a MDP and
solved by algorithms based on lookup tables (e.g., dynamic programming and
reinforcement learning). However, the state space required to capture the
channel and buffer status (see Sec. III-A) can be extremely large that
motivates us to develop a DRL-based approach.

##  III DDQN-PER Relay Selection Scheme

To minimize the average AoI, we propose a centralized model-free method based
on deep Q network (DQN), which leverages neural networks in place of look-up
tables to find $\boldsymbol{\mathcal{R}}$ while improving the algorithm
convergence and stability through double DQN (DDQN) with priority experience
replay (PER). We refer the proposed scheme as DDQN-PER.

###  III-A State Space

In making the RS decision, the agent at $S$ learns the environment from the
following four pieces of information: the CSI of the $S-R$ links and the $R-D$
links, the data and the energy buffer status. Firstly, the instantaneous CSI
is parameterized by the normalized channel gains as given by
$\bar{\boldsymbol{\mathcal{H}}}(n)=\frac{\boldsymbol{\mathcal{H}}(n)}{\max\;\boldsymbol{\mathcal{H}(n)}}$
where $\boldsymbol{\mathcal{H}}(n)=[|h_{1}(n)|~{}\cdots~{}|h_{R_{K}}(n)|]^{T}$
for the $S-R$ links and
$\bar{\boldsymbol{\mathcal{G}}}(n)=\frac{\boldsymbol{\mathcal{G}}(n)}{\max\;\boldsymbol{\mathcal{G}(n)}}$
where $\boldsymbol{\mathcal{G}}(n)=[|g_{1}(n)|~{}\cdots~{}|g_{R_{K}}(n)|]^{T}$
for the $R-D$ links. Provided with CSI, the agent can make a decision to meet
constraint (3d) while maximizing the total amount of harvested energy by
relays. For example, if two relays satisfy constraint (3d), the one with a
higher channel gain of the $S-R$ link could charge more energy and thus it can
be selected to replenish energy while the other relay is selected to receive a
packet from $S$. Clearly, the amount of accumulated energy depends on the
available room of the energy buffer. For slot $n$, the energy buffer status is
represented by
$\bar{\boldsymbol{\mathcal{E}}}(n)=\frac{\boldsymbol{\mathcal{E}}(n)}{E_{\max}}$
where
$\boldsymbol{\mathcal{E}}(n)=[\epsilon_{1}(n)~{}\cdots~{}\epsilon_{R_{K}}(n)]^{T}$.
Knowing $\bar{\boldsymbol{\mathcal{E}}}(n)$ is also necessary to meet the
energy buffer constraint in (3c). Similarly, the data buffer constraint in
(3c) is handled by feeding the data buffer status denoted as
$\boldsymbol{\mathcal{D}}(n)=[\delta_{1}(n)~{}\cdots~{}\delta_{R_{K}}(n)]^{T}$
to the agent. When more than one relays satisfy the constraints (3c) and (3d),
the age information is essential for decision-making. Therefore, we include
the normalized relative age denoted as
$\bar{\boldsymbol{\mathcal{A}}}(n)=\frac{\boldsymbol{\alpha^{\prime}}(n)}{A_{\max}}$
in the state, where
${\boldsymbol{\alpha}}^{\prime}(n)=[\alpha^{\prime}_{1}(n),\cdots,\alpha^{\prime}_{K}(n)]^{T}$
is the vector of relative age. In addition, a binary indicator denoted as
$I(n)$ is employed to indicate whether there is a packet arrives in slot $n$,
which takes the value of one when there is a new packet generated in slot $n$
and it is zero otherwise. Overall, the state of slot $n$ is represented as
$s(n)\triangleq\\{\bar{\boldsymbol{\mathcal{H}}}(n),\bar{\boldsymbol{\mathcal{G}}}(n),\boldsymbol{\mathcal{D}}(n),\bar{\boldsymbol{\mathcal{E}}}(n),\bar{\boldsymbol{\mathcal{A}}}(n),I(n)\\}$.
While $\boldsymbol{\mathcal{D}}(n)$, $\bar{\boldsymbol{\mathcal{E}}}(n)$ and
$\bar{\boldsymbol{\mathcal{A}}}(n)$ are finite in dimensions,
$\bar{\boldsymbol{\mathcal{H}}}(n)$ and $\bar{\boldsymbol{\mathcal{G}}}(n)$
have infinitely large dynamic ranges. Consequently, the dimension of $s(n)$ is
infinite.

###  III-B Action Space

At the beginning of each slot, the agent picks an action to determine which of
the $2N$ links is activated. Let $R^{*}(n)$ denote the relay associated with
the selected link in slot $n$, where $R^{*}(n)$ is a transmitting relay if an
$R-D$ link is activated or a receiving relay when an $S-R$ link is activated,
whichever maximizes the reward as will be specified later. Accordingly, the
action of slot $n$ can be expressed as $a(n)=R^{*}(n)$ and the action space of
$N$ slots is given by $\mathcal{A}=\\{a(1),\cdots,a(N)\\}$.

###  III-C Reward Function

Based on $s(n)$, the agent takes a specific action and receives a certain
reward according to a reward function. To minimize the AoI, it is better to
maximize packet freshness at the relay side. Noted that based on the
definition of relative age in Sec. II-A, the positive value of relative age
means the packet in the relay buffer is fresher than the packet received at
$D$ and the negative value of relative age means the packet is outdated. Thus,
the action of activating the $S-R$ link is assigned with the reward equal to
the accumulated relative age of all the packets in the relay buffers, as given
by

$\displaystyle r_{a}(n)=\sum_{k=1}^{K}\alpha_{k}^{{}^{\prime}}(n).$ (4)

As for the $R-D$ link, the action of transmitting the freshest packet should
be encouraged. This can be achieved by a reward function given as

$\displaystyle r_{a}(n)=\alpha^{\prime}_{R_{t}}(n).$ (5)

We note that some actions may be infeasible, when (i) any of the constraints
(3b)-(3e) are violated, (ii) there is a packet generated at $S$ but the action
of activating the $R-D$ link is chosen, and (iii) there is no packet generated
at $S$ but the action of activating the $S-R$ link is chosen. To remedy, an
infinitely small negative reward is assigned to the infeasible action.

###  III-D Proposed DDQN-PER-based RS Algorithm

In the proposed DDQN-PER-based RS algorithm, the agent employs two
convolutional neural networks (CNNs), namely the main network $Q$ and the
target network $\hat{Q}$ [15], for determining the predicted Q-value and the
target Q-value, respectively. These two networks share the same structure but
with different parameters. The main network updates its parameters denoted as
$\theta$ in each step while the target network synchronizes its parameters
denoted as $\theta^{-}$ with $\theta$ of the main network periodically. In the
training procedure, the agent executes the action determined by the maximum
Q-value of the main network that is updated in every step to explore the
unknown environment. Specifically, the system state changes from the current
state $s(n)$ to the next state $s(n+1)$ when the agent selects an action
$a(n)$ and earns an reward $r(n)$. The target Q-value in slot $n$ is

$\displaystyle
y(n)=r(n+1)+\omega~{}[\hat{Q}(s(n+1),\arg\max_{a(n)}Q(s(n+1),a(n))]$ (6)

where $\omega\in[0,1]$ is the discount factor and $Q(\cdot)$ and
$\hat{Q}(\cdot)$ represent the action-value function of main network and
target network, respectively. The mean square error between the target value
and the predicted value of all transitions is defined as the loss function
given by

$\displaystyle
L(\theta)=E\bigl{[}y(n)-Q\bigl{(}s(n),a(n)\bigr{)}\bigr{]}^{2}.$ (7)

By performing gradient descent to $L(\theta)$, the main network updates its
parameters in each step.

###  III-E Prioritized Experience Replay

Instead of sampling the experience uniformly, more efficient learning can be
achieved by sampling the experience according to its importance, which is
measured by the magnitude of the loss function of each transition. Let
$L_{i}(\theta)$ denote the loss function of transition $i$, which is assigned
with priority $p_{i}=|L_{i}(\theta)|+\zeta$, with $\zeta$ being a small
positive constant to avoid the transition with a zero loss function value
never getting sampled. Given $p_{i}$, the probability of transition $i$ is

$\displaystyle P_{i}=\frac{p_{i}^{\alpha}}{\sum_{j}p_{j}^{\alpha}}$ (8)

where $\alpha\geq 0$ is used to control the weight of priority in the sampling
process. The case $\alpha=0$ corresponds to uniform sampling. Clearly,
prioritized replay changes the sample distribution that leads to biased
estimate of the solution. The bias can be fixed by applying importance-
sampling weights as given by [16]

$\displaystyle w_{i}=(\frac{1}{\mathcal{RM}}\cdot\frac{1}{P_{i}})^{\beta}$ (9)

where $\mathcal{RM}$ is the replay memory size and $\beta$ is an exponent that
controls how much prioritization to apply. The case $\beta=1$ reduces to
uniform sampling without prioritization.

The complete DDQN-PER RS algorithm is presented in Algorithm. 1. In Line 6-8,
action $a(n)$ is selected based on $\varepsilon$-greedy and the corresponding
experience is stored in the replay memory. Line 9-15 calculate the loss
function, the importance-sampling weight, and transition probability for each
transition in the minibatch and accumulate the loss function value with weight
$\theta$. In Line 16, the main network is updated. Lastly, the target network
is periodically updated every $C$ steps.

Algorithm 1 DDQN-PER-based relay selection algorithm

1: Input: mini-batch size $\mathcal{M}$, step-size $\mu$, replay memory size
$\mathcal{RM}$, exponents $\alpha$ and $\beta$.

2: Initialize replay memory, $\Delta=0$, $p_{1}=1$

3: for $\text{Episode}=1$ to $E$ do

4: Initialize each data buffer and energy buffer.

5: for $n=1$ to $N$ do

6: Select $a(n)=\arg\max_{a}\mathcal{Q}\allowbreak(s(n),a)$ with probability
$\varepsilon$ or a random action with probability $1-\varepsilon$.

7: Execute action $a(n)$ and observe reward $r(n)$ and the corresponding state
$s(n+1)$

8: Store $z_{n}=\\{s(n),a(n),r(n),s(n+1)\\}$ with maximal priority
$p_{n}=\max_{j<n}p_{j}$ in replay memory

9: for $i$=1 to $\mathcal{M}$ do

10: Sample transition $i$ by probability $P_{i}$ in (8).

11: Compute $w_{i}$ using (9).

12: Compute loss function value $L_{i}(\theta)$ in (7)

13: Update transition priority $p_{i}\leftarrow|\rho_{i}|+\zeta$

14: Accumulate weight change
$\Delta\leftarrow\Delta+w_{i}L_{i}(\theta)\Delta_{\theta}\times
Q(s_{i},a_{i})$

15: end for

16: Update weights $\theta\leftarrow\theta+\mu\cdot\Delta$, reset $\Delta=0$

17: Reset $\theta^{-}=\theta$ every $C$ steps

18: end for

19: end for

##  IV Results and Discussions

###  IV-A Simulation Setting

In implementing DDQN, the main network $Q$ and the target network $\hat{Q}$
have the same structure comprising of an input layer with dimension $K\times
4$, a fully connected layer with eighty nodes using a rectified linear unit
(ReLU) activation function, and an output layer with the dimension of
$|\mathcal{A}|\times 1$. The DQN is trained using 1,000 episodes with 2,000
steps in each episode. The discount factor $\omega$ is set to one, the
prioritized factors $\alpha=0.6$ and $\beta=0.4$. The size of the replay
memory $\mathcal{RM}$ and mini-batch $\mathcal{M}$ is 10,000 and 32,
respectively. The term “SNR” in the figures is defined as
$\text{SNR}=P_{S}/N_{0}$, which is set to 70 dB, the noise power $N_{0}=0.01$,
and transmission power ratio $P_{S}/P_{R}=1,000$. For correct decoding, the
SNR threshold is equal to $2^{2R-1}$ where $R=1$ bit/sec/Hz is the target
transmission rate. The packet arrival probability is $\lambda=0.3$ and
$A_{\max}=100$. $E_{\max}=3$ and the RF-DC conversion efficiency $\eta=0.5$
for all the relays. Initially, each data buffer is empty, and each energy
buffer is filled with one energy interval. All the links have the same
distance of 36 m.

###  IV-B Benchmarks

We compare the performance of the proposed DQN-PER RS scheme with three
benchmark schemes.

  * •

Dual Buffer Relay Selection (DBRS) [12], which is a heuristic scheme for
minimizing the outage probability of buffer-aided relay networks. This scheme
is chosen to verify whether minimizing the outage probability is beneficial to
AoI minimization and it is similar to the proposed scheme with both CSI and
buffer status considered.

  * •

Greedy scheme, which makes the RS decision based on the generation time of
each packet. The relay that has the latest generation time is selected to
transmit. As to the reception relay, the relay with an empty data buffer has
the highest priority to receive a status packet from $S$. If more than one
relay has an empty buffer, the tight is broken by a random selection. If more
than one relay has a non-empty buffer, the relay with the earliest generated
packet is selected as $R_{r}$.

  * •

Max-Link scheme [17], which activates a link that has the highest SNR among
all eligible $S-R$ and $R-D$ links, where an eligible link must satisfy both
constraints (3c) and (3d). This selection rule is similar to the proposed
scheme but none of the buffer status is considered.

We note that DBRS and Greedy activate the $S-R$ link and the $R-D$ link in a
fixed order while the order can be arbitrary in Max-Link and the proposed
schemes. The activation order would affect the status update performance as
will be shown later. Also, if the selected relay by DBRS and Greedy schemes
violates constraints (3c) or (3d), a single relay is chosen as:
$R^{*}=\max_{\epsilon_{R_{k}}\neq
0,\delta_{R_{k}}<D_{\max}}~{}\min\\{\gamma_{h_{k}},\gamma_{g_{k}}\\}$ [12].

###  IV-C Performance evaluation

Fig. 1 shows the training result of the proposed DDQN-PER model by plotting
the loss function values for 500 episodes. One can clearly observe that the
proposed DRL algorithm converges for varied number of relays $K$. For a
greater $K$, the loss function curve tends to converge ealrier.

Figure 1: Training loss of the proposed DDQN-PER model.

Fig. 2 depicts the average AoI of different RS schemes versus the number of
relays $K$. As $K$ increases, the average AoI of all the RS schemes is
effectively reduced, thanks to the increased selection diversity. It is
noticed that the proposed scheme has a higher packet lost rate (the ratio
between the number of discarded packets by relays and the total number of
packets generated by $S$) than _DBRS_ and _Greedy_ schemes but its average AoI
is significantly lower. However, the same gain is not attained using _Max-
Link_ scheme, which has a higher packet loss rate. This suggests that properly
dropping packets helps to improve the AoI performance. For all the schemes,
the gain from increasing $K$ is minor. This can be explained because the
status packet is backlogged at relays when none of the relays has a sufficient
amount of energy to transmit and this chance does not diminish as $K$
increases [12]. It is worth to point out _Max-Link_ scheme performs the worst
when $K$ is large since it chooses to activate a link without considering
packet freshness. _DBRS_ also does not consider packet freshness but it
performs better than _Max-Link_ by further considering buffer status. Although
_Greedy_ scheme considers packet freshness, it lacks of considering channel
and buffer status that might end in choosing an infeasible action.

Figure 2: Average AoI and packet loss rate vs. number of relays $K$.

The impact of energy buffer size $E_{\max}$ is investigated in Fig. 3 for
$K=3$ relays. Intuitively, a larger energy buffer helps to accumulate more
energy thereby less likely that none of relays can transmit. From the figure,
only marginal reduction of the average AoI is achieved by increasing the
energy buffer size since the amount of harvested energy is relatively small
compared to that required for transmitting a packet. In other words, the
energy consumption rate is much higher than the supply rate and this unbalance
can not be mitigated by a large energy buffer size.

Figure 3: Average AoI vs. energy buffer size $E_{\max}$ for $K=3$ relays.

##  V Conclusion

This work addressed the sophisticated RS problem for AoI minimization in the
EH relay-assisted status update network by proposing a DQN based approach.
Simulation results demonstrated that the proposed RS scheme effectively
improves the average AoI than the other competitive schemes. When the number
of relays is small, the proposed scheme reduces about $50\%$ average AoI of
existing RS schemes. Among numerous factors that affect the AoI performance,
adding many relays and the size of energy buffers only has limited improvement
to the average AoI.

## References

  * [1] W.-Y. Chung and et al., "A cell phone based health monitoring system with self analysis processor using wireless sensor network technology," _Proc. 29th Annual International Conference of the IEEE Engineering in Medicine and Biology Society_ , 2007, pp. 3705-3708.
  * [2] H. Liu, Z. Meng, and S. Cui, “A wireless sensor network prototype for environmental monitoring in greenhouses,” in _Proc. International Conference on Wireless Communications, Networking and Mobile Computing_ , 2007, pp. 2344–2347.
  * [3] 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. 2, pp. 1807–1827, 2019.
  * [4] M. Moltafet, M. Leinonen, and M. Codreanu, “On the age of information in multi-source queueing models,” _IEEE Trans. Commun._ , vol. 68, no. 8, pp. 5003–5017, 2020.
  * [5] L. Liu, K. Xiong, J. Cao, Y. Lu, P. Fan, and K. B. Letaief, “Average AoI minimization in UAV-assisted data collection with RF wireless power transfer: A deep reinforcement learning scheme,” _IEEE Internet Things J._ , vol. 9, no. 7, pp. 5216–5228, 2022.
  * [6] A. Bletsas, A. Khisti, D. Reed, and A. Lippman, “A simple cooperative diversity method based on network path selection,” _IEEE J. Sel. Areas Commun._ , vol. 24, no. 3, pp. 659–672, 2006.
  * [7] M. Moradian and A. Dadlani, “Age of information in scheduled wireless relay networks,” in _Proc. 2020 IEEE Wireless Communications and Networking Conference (WCNC)_ , 2020, pp. 1–6.
  * [8] M. Xie, J. Gong, and X. Ma, “Age-energy tradeoff in dual-hop status update systems with the m-th best relay selection,” in _Proc. IEEE Vehicular Technology Conference (VTC2021-Spring)_ , 2021, pp. 1–5.
  * [9] ——, “Age and energy tradeoff for short packet based two-hop decode-and-forward relaying networks,” in _Proc. IEEE Wireless Communications and Networking Conference (WCNC)_ , 2021, pp. 1–6.
  * [10] T. D. P. Perera, D. N. K. Jayakody, I. Pitas, and S. Garg, “Age of information in SWIPT-enabled wireless communication system for 5GB,” _IEEE Wireless Commun._ , vol. 27, no. 5, pp. 162–167, 2020.
  * [11] D. Khan and et al., “A survey on RF energy harvesting system with high efficiency RF-DC converters,” _Journal of Semiconductor Engineering_ , vol. 1, no. 1, pp. 13–30, 2020.
  * [12] C.-H. Lin and K.-H. Liu, “Relay selection for energy-harvesting relays with finite data buffer and energy storage,” _IEEE Internet Things J._ , vol. 8, no. 14, pp. 11 249–11 259, 2021.
  * [13] A. Zakeri, M. Moltafet, M. Leinonen, and M. Codreanu, “Minimizing AoI in resource-constrained multi-source relaying systems with stochastic arrivals,” in _Proc. IEEE Global Communications Conference (GLOBECOM)_ , 2021, pp. 1–6.
  * [14] R. D. Yates, Y. Sun, D. R. Brown, S. K. Kaul, E. Modiano, and S. Ulukus, “Age of information: An introduction and survey,” _IEEE J. Sel. Areas Commun._ , vol. 39, no. 5, pp. 1183–1210, 2021.
  * [15] H. Van Hasselt, A. Guez, and D. Silver, “Deep reinforcement learning with double Q-learning,” in _Proc. AAAI Conference on Artificial Intelligence_ , vol. 30, no. 1, 2016.
  * [16] T. Schaul, J. Quan, I. Antonoglou, and D. Silver, “Prioritized experience replay,” _arXiv:1511.05952_ , 2015.
  * [17] I. Krikidis, T. Charalambous, and J. S. Thompson, “Buffer-aided relay selection for cooperative diversity systems without delay constraints,” _IEEE Trans. on Wireless Commun._ , vol. 11, no. 5, pp. 1957–1967, 2012.


# Re-understanding of the deformation potential constant in the single crystal
silicon

Aijun Hong<EMAIL_ADDRESS><EMAIL_ADDRESS>School of Physics,
Communication and Electronics, Jiangxi Normal University, Nanchang 330022,
China  Feng Sun  School of Education, Nanchang Institute of Science and
Technology, Nanchang 330108, China

###### Abstract

The mobility formula based on deformation potential (DP) theory is of great
importance in semiconductor physics. However, the related calculations for the
DP constant are controversial. It is necessary to redo in-depth and
comprehensive research on the mobility of single crystal silicon and the
related parameters such as the effective mass and the DP constant. In this
work the mobility effective mass is defined and a method based on the first
principles is presented to evaluate the correction of the DP constant. It is
found that the effective mass is closer to experimental data and the
correction of the DP is a negligible value of about 0.3 eV being far lower
than $\sim$10 eV calculated by other methods. Using these parameters, we
obtain the mobilities of the single crystal silicon in reasonable agreement
with the experimental values.

##  I Introduction

Since the deformation potential (DP) theory was proposed by Bardeen and
Shockley in 1950 [1], it has been applied to describe the interaction between
electron and phonon [2, 3, 4] and evaluate transport properties of
semiconductors successfully [5, 6, 7]. In particular, its combination with
density functional theory (DFT) provides a favorable mean for explaining
experiments and predicting new materials in the thermoelectric community.
However, there are still some disagreements about the calculation for the DP
constant of bulk materials. The focus of disagreements is that no precise
reference energy level is found for obtaining the absolute band energy of the
different deformation structures.

As is well known, the effective one-electron Schrödinger (or Dirac) equation
with respect to the eigenstates and eigenvalues ($\Psi_{i}$ and
$\varepsilon_{i}$) is expressed as [8]

$[K_{op}+V_{eff}(\textbf{r})]\psi_{i}(\textbf{r})=\epsilon_{i}\psi_{i}(\textbf{r}),$
(1)

where the effective Coulomb potential operator is the sum of the exchange
correction potential ($\mu_{xc}$) and Coulomb potential ($V_{c}$)

$V_{eff}(\textbf{r})=\mu_{xc}(\textbf{r})+V_{c}(\textbf{r}).$ (2)

The $V_{c}$ is given by

$V_{c}(\textbf{r})=\int\frac{\rho(\textbf{r}^{\prime})}{|\textbf{r}-\textbf{r}^{\prime}|}d\textbf{r}^{\prime}-\sum_{\alpha}\frac{Z_{\alpha}}{|\textbf{r}-\textbf{R}_{\alpha}|}.$
(3)

Obviously, there are the Coulomb singularities at the electron and nucleus
sites ($\textbf{r}^{\prime}$ and $\textbf{R}_{\alpha}$). In order to solve the
Eq. (1), there are many methods proposed to deal with the singularities.
However, it is impossible to avoid introducing an uncertain constant when
solving the Poisson’s equation in any case. In fact, introducing the uncertain
constant ignores the effects of both volume deformation and atomic species on
the periodic potential. As a result, it makes the band energies of different
materials and of the same material with different volumes, unable to be
compared quantitatively, although the influence on the total energy can be
canceled by special mathematical processing [9, 10, 11].

For the calculation of the DP constant of one material, it is not necessary to
consider the influence of atomic species. Therefore, finding a physical
quantity that is not affected by the volume deformation as the reference
energy level is the key to calculating the DP constant. At present, there are
average potential [12, 13, 14, 15] and core state level [16, 17] alignments to
define the reference energy level. The former is also known as the supercell
method, where the band structures are aligned according to the average
potential from the core electrons in the supercell composed of the compressed
and expanded parts. The latter is based on the assumption that the energy
level of the core electron is hardly affected by the volume deformation.

Usually, the core state level alignment method is adopted to explore
thermoelectric properties due to the requirement of less computational
resources. In fact, the core state levels are sensitive to the volume
deformation, and thus have non-negligible DP values [16]. Consequently, the DP
constants corrected by the core state level alignment are small and thus lead
to much high mobility deviating from the experimental results. The
disagreement between theory and experiment cannot be completely attributed to
the theoretical calculation, because the recognized accurate experimental
method for determining the DP constant has not been proposed.

Experimentally, the present determinations of the DP constant are only
performed indirectly with a significant amount of assumptions. Therefore, the
experimental DP constant for the same material has great differences. Taking
GaAs as an example, Nolte et al. used the impurity level as the reference
level to obtain the DP value of -9.3 eV for the conduction band minimum (CBM)
[18], however, its value is 7.0 eV from the mobility measurements at low
temperature [19]. Interestingly, a value of $16.0\pm 0.5$ eV for CBM is
obtained from high-temperature mobility [20]. Moreover, Pfeffer et al. [21]
inferred from the optical absorption in Se-doped GaAs bulk [22] that the DP
value is as high as 15.7 eV. A much higher value of 20 eV is accepted and
adopted in the Supporting Information of the recent literature [23]. In a
word, the experiment has not so far given satisfactory results about the DP
constant for GaAs.

Theoretically, the DP constant of a material also varies greatly. In the
original literature of DP theory [1], Bardeen and Shockley first deduced 11.3
eV for the DP value of the valence band maximum (VBM) of cubic silicon from
measured mobility data in the polycrystalline silicon. Latter, the values
-10.2 eV [24], -7.9 eV [25], and 2.05 eV [17] are proposed by different
researchers. It was claimed in some literatures [14, 16, 26, 27] that the
obtained $\sim$1.5 eV along the [100] direction is very consistent with the
experimental value $1.8\pm 0.7$ eV [28] for the p-type silicon. However, it is
puzzled that we have not found relevant experimental data and only see the
value of DP constants ($3.3\pm 0.7$ eV) for the n-type silicon in Ref. [28].
Furthermore, it is worth mentioning that the experimental data is obtained in
the heavily doped n-type silicon with the high concentration of $5\times
10^{21}cm^{-3}$ and its Fermi level is located above the CBM. Therefore, the
measured DP constant is for the Fermi level above the CBM, however, the
theoretical DP constant is for the CBM [14]. This implies that the theory and
experiment are not consistent, even may deviate.

In the present work we first redefined the relevant physical parameters in the
mobility formula, based on the analysis and summary of previous theoretical
and experimental data on the DP constant. Then we present a strategy to
evaluate the correction of the DP constants of materials with the aid of all-
electron first-principles approach. Finally, we employ the modified mobility
formula to calculate the mobility of the single crystal silicon for verifying
the effectiveness of the strategy.

##  II Methodology

According to the DP theory, the carrier mobility $\mu$ of non-polar crystals
with the properties of isotropy and non-degeneracy is expressed as [1]

$\mu_{e(h)}=\frac{2\sqrt{2\pi}e\hbar^{4}\rho
v^{2}}{3(k_{B}T)^{3/2}m^{5/2}_{e(h)}(\lambda_{e(h)})^{2}},$ (4)

where $k_{B}$, $e$ and $\hbar$ represent, respectively, the Boltzmann
constant, the electronic charge and the reduced Planck constant, and
$\rho$,$v$, $m_{e(h)}$ and $\lambda_{e(h)}$ are the mass density, the
longitude acoustic velocity, the effective mass and DP constant for the n- or
p-type (electron/hole) carrier (namely, for the CBM or valence band maximum
(VBM)).

The $v$ of the isotropic polycrystal is govern by [29]

${{\textit{v}}}={{\left(\frac{3{{B}}+4G}{3\rho}\right)}^{1/2}},$ (5)

where $B$ and $G$ are the bulk modulus and the shear modulus, respectively.
Sometimes the following expression, neglecting the shear strain contribution
to the longitude acoustic velocity, is adopted.

${{\textit{v}}}={{\left(\frac{{{B}}}{\rho}\right)}^{1/2}}.$ (6)

This is reasonable because small shear strain does not induce the volume
change and therefore cannot generate effective DP value.

The longitude acoustic velocity $v^{\beta}$ along the $\beta$ direction in the
single crystal is approximatively expressed by the elastic constant
$c_{\beta}$ and the mass density,

$v^{\beta}=\sqrt{\frac{c^{\beta}}{\rho}}.$ (7)

For the anisotropic single crystal, the effective mass of electron/hole along
the $\beta$ direction ($m^{\beta}_{e(h)}$) can be obtained by calculating the
band effective mass at the CBM/VBM

$\frac{1}{m^{\beta}_{e(h)}}=\frac{\partial^{2}\varepsilon}{\hbar^{2}\partial(k^{\beta})^{2}}|_{\varepsilon=\varepsilon_{CBM(VBM)}},$
(8)

where $\epsilon$ and $k^{\beta}$ stand for the energy of band edge and the
wave vector along the $\beta$ direction.

We define three types of DP constants

$\lambda_{e(h)}^{X}=\frac{\mathrm{d}\varepsilon}{\mathrm{d}\delta^{X}}|_{\varepsilon=\varepsilon_{CBM(VBM)}},$
(9)

where $\delta^{X}$ represents different types of strains ($X=L,S\ or\ V$ for
the uniaxial, biaxial or volumetric strain), $\varepsilon_{CBM(VBM)}$ denotes
the energy at the CBM/VBM

$\displaystyle\delta^{L}=\frac{L-L_{0}}{L_{0}},$ (10)
$\displaystyle\delta^{S}=\frac{S-S_{0}}{S_{0}},$ (11)
$\displaystyle\delta^{V}=\frac{V-V_{0}}{V_{0}},$ (12)

where $L_{0}$, $S_{0}$ and $V_{0}$ ($L$, $S$ and $V$) are length, area and
volume of crystals without (with) the strains. Note that the DP constant
$\lambda^{\beta}_{e(h)}$ can be obtained by applying the uniaxial strain, and
thus is equivalent to the $\lambda^{L}_{e(h)}$. Also, we call
$\lambda^{X}_{e(h)}$ ($L$, $S$ and $V$) as the uniaxial, biaxial and
volumetric DP constant, respectively.

Usually, the CBM/VBM of the crystal with high symmetry has many equivalent
points that possess the same energy and commonly determine the transport
properties. The points can be divided into two categories (copoint or non-
copoint degeneracy) based on whether they are at the same position. For
instance, the n-type single crystal silicon has six non-copoint degenerate
CBMs which symmetrically locate at the _a_ *, _b_ * and _c_ * axes of the
reciprocal space. Each two CBMs on the same axis have the same electronic
structural properties along the same direction.

Each CBM has an ellipsoidal isoenergetic surface near it, and independently
participates in electrical transport. Moreover, each CBM exhibits anisotropic
transport behavior because the effective mass varies along different
directions of the ellipsoidal isoenergetic surface. However, the total
transport of the six CBMs is isotropic. For _i_ th CBM, we assume that its
contribution to electrical conductivity along one special $\beta$ direction
($\sigma_{i}^{\beta}$) is equivalent to the contribution of a spherical
isoenergetic surface with the effective mass of the ellipsoidal isoenergetic
surface along the $\beta$ direction, rather than the assumption of constant
relaxation time. Thus, the total electronic conductivity along the $\beta$
direction can be expressed as

$\sigma^{\beta}=\sum_{i=1}^{N_{v}}{\sigma^{\beta}_{i}}=\sum_{i=1}^{N_{v}}{n^{\beta}_{i}qu^{\beta}_{i}}={nqu^{\beta}},$
(13)

where $N_{V}$, $u^{\beta}$ and _n_ are the number of band degeneracies at band
edge, the total real mobility and the total real carrier density that is
described by

$n=\sum_{i=1}^{N_{v}}{n_{i}}=\frac{V_{c}}{2\pi^{2}}(\frac{2m_{i}}{\hbar^{2}})^{\frac{3}{2}}\varepsilon^{\frac{1}{2}}f,$
(14)

Here, $n_{i}$ is the real carrier density supported by the _i_ th band
degeneracy, $\varepsilon$, $V_{c}$ and $f$ are the energy of band edge, the
volume of unit cell and the Fermi-Dirac distribution function, and $m_{i}$ is
the local DOS effective mass of _i_ th band degeneracy. The $m_{i}$ is
determined by the shape of the isoenergetic surface near the band edge. The
$m_{i}$ of _i_ th spherical isoenergetic surface is band effective mass at
CBM/VBM, which is equal to the effective mass along arbitrary $\beta$
direction ($m_{i}^{\beta}$). For elliptical isoenergetic surface, it can be
written as

$m_{i}=\sqrt[3]{m_{i}^{x}m_{i}^{y}m_{i}^{z}},$ (15)

where $m_{i}^{x},m_{i}^{y}$ and $m_{i}^{z}$ are the effective masses along the
three principal axes of the ellipsoidal isoenergetic surface. The $m_{i}$ is
the result of equating the volume of an ellipsoid to a sphere. Mathematically
speaking, it is the geometric mean of the effective masses along the three
direction. The $m_{i}$ expression for the case of anisotropy has significant
limitations. For example, when the effective mass along the one direction is
infinitely close to zero, the $m_{i}$ is close to zero. Obviously, this is not
in line with the facts. Therefore, it seems more reasonable to use geometric
mean of the effective masses to describe the DOS near band edge.

$m_{i}=\frac{m_{i}^{x}+m_{i}^{y}+m_{i}^{z}}{3},$ (16)

In Eqs.(13), $n^{\beta}_{i}$ and $u^{\beta}_{i}$ are the equivalent carrier
density and the equivalent mobility of _i_ th band degeneracy given by

$\displaystyle n^{\beta}_{i}$
$\displaystyle=\frac{V_{c}}{2\pi^{2}}(\frac{2m^{\beta}_{i}}{\hbar^{2}})^{\frac{3}{2}}\varepsilon^{\frac{1}{2}}f,$
(17)

$\displaystyle\mu_{i}^{\beta}$
$\displaystyle=\frac{2\sqrt{2\pi}e\hbar^{4}c^{\beta}}{3(k_{B}T)^{3/2}(m_{i}^{\beta})^{5/2}(\lambda^{\beta})^{2}}.$
(18)

Here, $m^{\beta}_{i}$ represents the effective masses of _i_ th band CBM/VBM
along the same $\beta$ direction. Substituting Eqs. (14), (17) and (18) into
(13), the total real mobility is obtained

$\displaystyle\mu^{\beta}$
$\displaystyle=\frac{2\sqrt{2\pi}e\hbar^{4}c^{\beta}}{3(k_{B}T)^{3/2}(m_{[m]}^{\beta})^{5/2}(\lambda^{\beta})^{2}}.$
(19)

We define the $m_{[m]}^{\beta}$ as the mobility effective mass along the
$\beta$ direction,

$\displaystyle\frac{1}{(m_{[m]}^{\beta})^{5/2}}=\frac{\sum\limits_{i=1}^{N_{v}}\frac{1}{m_{i}^{\beta}}}{\sum\limits_{i=1}^{N_{v}}(m_{i})^{\frac{3}{2}}}.$
(20)

When the $N_{v}$ band degeneracies have the same isoenergetic surfaces, one
can get the more common variation of Eq. (20),

$\displaystyle\frac{1}{(m_{[m]}^{\beta})^{5/2}}=\frac{\frac{1}{N_{v}}\sum\limits_{i=1}^{N_{v}}\frac{1}{m_{i}^{\beta}}}{\frac{1}{N_{v}}\sum\limits_{i=1}^{N_{v}}(m_{i})^{\frac{3}{2}}}=\frac{\frac{1}{m_{c}}}{m_{i}^{3/2}}=\frac{1}{m_{c}m_{i}^{3/2}}.$
(21)

The $m_{c}$ can be called the conductivity effective mass. In fact, it is the
mathematical harmonic mean. For the weak anisotropy, the ratio of $m_{c}$ and
$m_{i}$ is close to 1. Thus, the other vibration[30] of Eq. (20) is

$\displaystyle\frac{1}{(m_{[m]}^{\beta})^{5/2}}\approx\frac{1}{m_{i}^{5/2}}.$
(22)

Moreover, when the same isoenergetic surfaces locate at the same position,
intervalley scattering can proceed smoothly like intravalley scattering.
Therefore, the Eq. (18) should be writen as

$\displaystyle\mu_{i}^{\beta}$
$\displaystyle=\frac{2\sqrt{2\pi}e\hbar^{4}c^{\beta}}{3N_{V}(k_{B}T)^{3/2}(m_{i}^{\beta})^{5/2}(\lambda^{\beta})^{2}}.$
(23)

and the Eq. (20) has a variant formula

$\displaystyle\frac{1}{(m_{[m]}^{\beta})^{5/2}}=\frac{\sum\limits_{i=1}^{N_{v}}\frac{1}{m_{i}^{\beta}}}{N_{V}\sum\limits_{i=1}^{N_{v}}(m_{i})^{\frac{3}{2}}}=\frac{1}{N_{V}m_{c}m_{i}^{3/2}}$
(24)

For n-type silicon, when the $\beta$ direction is parallel to the $a^{*}$
axis, the $m^{\beta}_{1}$=$m^{\beta}_{4}$, $m^{\beta}_{2}$=$m^{\beta}_{5}$ and
$m^{\beta}_{3}$=$m^{\beta}_{6}$ correspond to $m_{l}$, $m_{t}$ and $m_{t}$
calling the longitude and transverse effective masses. Six $m_{i}$ values is
equal and it is given by

$\displaystyle m_{i}=\sqrt[3]{m_{t}m_{t}m_{l}}.$ (25)

or

$\displaystyle m_{i}=\frac{m_{t}+m_{t}+m_{l}}{3}.$ (26)

Note that the effect of the intervalley scattering on the mobility is not
taken account into in n-type silicon, it is because the CBMs possess vastly
different momentum (namely locate at different sites), which hinders the
intervalley scattering.

Different from the n-type silicon, the p-type silicon has three VBMs
converging at $\Gamma$ point. The three degenerate bands contain two heavy
bands and one light bands, so $m_{i}^{\beta}$ (i=1, 2, 3) is labled by
$m_{hh}$, $m_{hh}$ and $m_{lh}$ respectively calling heavy and light band
effective masses. Moreover, the carrier transition occurs between the same
momentums and between the same energies, namely, the scattering near the VBMs
is approximatively elastic. We assumed that the intervalley scattering is as
well as the intravalley scattering. Therefore, the scattering rate is three
times the case without intervalley scattering, and the hole mobility of the
single crystal is given by

$\mu_{h}^{\beta}=\frac{2\sqrt{2\pi}e\hbar^{4}c^{\beta}}{3(k_{B}T)^{3/2}(m_{[m]h}^{\beta})^{5/2}(\lambda^{\beta}_{h})^{2}}.$
(27)

With the help of Eq. (28), the $m_{[m]h}^{\beta}$ is roughly expressed as

$\displaystyle\frac{1}{(m_{[m]h}^{\beta})^{5/2}}=\frac{\sum\limits_{i=1}^{3}\frac{1}{m_{i}^{\beta}}}{3\sum\limits_{i=1}^{3}(m_{i})^{\frac{3}{2}}}=\frac{1}{3m_{c}m_{i}^{3/2}}.$
(28)

In fact, its accurate expression is

$\displaystyle\frac{1}{(m_{[m]h}^{\beta})^{5/2}}=\frac{\sum\limits_{i=1}^{3}(m_{i}^{\beta})^{1/2}}{(\sum\limits_{i=1}^{3}(m_{i})^{3/2})^{2}}.$
(29)

It is worth mentioning that the $\lambda^{\beta}$ is replaced by the
$\lambda^{V}$ in the fact calculations because it is difficult to accurately
calculate the deformation potentials of different degenerate band edges.
Therefore, the anisotropy of the mobility comes from the contribution of the
anisotropy of effective mass and elastic constant that respectively
characterize the properties of electrons and phonons.

##  III Results and discussion

Figure 1: (Color online) The uniaxial/biaxial/volumetric strain
($\varepsilon_{L}$, $\varepsilon_{S}$ and $\varepsilon_{V}$) effects on the
total energy (a)-(c), band gap (d)-(f), VBM/CBM energy (g)-(i), and average
energy of the DV state (j)-(l).

The relative total energies $\Delta E$ of the strained silicon with reference
to the unstrained structure are all positive values shown in Figs. 1(a)–1(c).
This indicates that the relative error for calculated lattice constant should
be less than 1$\%$. It is expected that the calculated lattice constant 5.46 Å
for Si agrees well with the experimental value 5.42 Å [31].

The improvements of total energies for the strained structures are related to
the breaking of symmetry and the interaction between atoms. Compared to
volumetric strain, the uniaxial and biaxial strains have greater impacts on
total energy of the system, because the both lead to the reduction of
structural symmetry. The uniaxially and biaxially strained structures belong
to space group ($I4_{1}$, No.141), nevertheless, the volumetrically strained
structure maintains the same structural symmetry ($Fd\overline{3}m$, No.227)
with the unstrained structure. At the same strain strength, the compression
more than the tensility can affect the total energy of silicon, because the
repulsive interactions between atoms are more sensitive to the relative
positions than the attractive interactions. For instance, the $\Delta E$ at
the $-1\%$ uniaxial strain is 0.0025 eV obviously larger than 0.0012 eV at the
$1\%$.

The calculated bandgap of unstrained silicon is 1.229 eV, close to the
experimental value of 1.17 eV at 0 K [32]. The quit small bandgap deviating is
mainly due to the calculated lattice constant slightly larger than the
experimental data. The good agreement is because the band structure in our
work is obtained within the modified Becke-Johnson exchange potential and the
generalized gradient approximation for the correlation (mBJ-GGA) [33], using
the full potential method as implemented in the WIEN2k code [34]. However, the
mBJ-GGA strategy has a shortcoming that no exchange and correlation energy
functional is defined [35], so that it cannot be used for all calculations in
connection with the total energy of the system, such as total energy and
geometric optimization calculations. As a consequence, the total energies in
this work are obtained by the standard GGA strategy.

Figure 2: (Color online) (a) Temperature dependence of the Fermi level
$\varepsilon_{F}$ at fixed n-type doping density(5$\times$$10^{21}$$cm^{-3}$).
(b),(c) Dependence of the Fermi level $\varepsilon_{F}$ on the p- and n-type
doping density at 300 K.

It is found in Figs. 1(d) and 1(e) that under uniaxial/biaxial strain the
bandgap is positively related to the compressive volume, and negatively to the
expansive volume. It is because both the compressive and tensile strains have
different effects on the conduction and valence bands. However, the
compressive and tensile strains have almost similar effects on the band edge
so that the bandgap and volumetric strain are positively correlated throughout
the overall region (see Fig. 1(f)). Without any energy corrections in Figs.
1(g)–1(i), we found the energy of the band edge goes down as the volume
increases, and the standard linear relationship between the band edge energy
and the volumetric tensile/compressive strain. However, the CBM energy is more
sensitive to the uniaxial tensile strain, on the contrary, the VBM energy is
more sensitive to the uniaxial compressive strain, which is also the case in
the biaxial strain. The corresponding DP constants are summarized in Table 1.
Under compressive strain, the $|\lambda_{VBM}^{L}|$ is as high as 14.90 eV
that is about seven times of $|\lambda_{CBM}^{L}|$, and is twice
$|\lambda_{VBM}^{L}|$ under tensile strain. Moreover, the overall
$|\lambda_{VBM}^{X}|$ is larger than the overall $|\lambda_{CBM}^{X}|$.
Especially, the difference between the overall $|\lambda_{VBM}^{L}|$ and
$|\lambda_{CBM}^{L}|$ reaches as large as 4.59 eV.

Table 1: Summary of the DP constans (in units of eV) for VBM and CBM. The
superscripts $L$, $S$ and $V$ represent obtaining DP constants through
uniaxial, biaxial and volumetric strains, respectively. "Compressed" and
"Expanded" represent the DP constant obtained through linear fitting of the
data from negative strain, positive strain and the both. "1S" and "2S"
represent the corrected DP constant with 1S or 2S core state as the reference
energy level, respectively.

| Compressed | Expanded | This work | 1S | 2S | Other work
---|---|---|---|---|---|---
$\lambda^{L}_{VBM}$ | -14.90 | -6.91 | -10.92 | 1.53 | 0.66 |  VBM: -10.2111Taken from Ref. [24]., -7.9222Taken from Ref. [25], 2.05333Taken from Ref. [17] CBM: 3.1444Taken from Ref. [14], -5.6555Taken from Ref. [25]
$\lambda^{L}_{CBM}$ | -2.09 | -10.52 | -6.33 | 6.11 | 5.25
$\lambda^{S}_{VBM}$ | -10.96 | -6.80 | -8.85 | 3.60 | 2.73
$\lambda^{S}_{CBM}$ | -6.27 | -10.38 | -8.31 | 4.14 | 3.27 |
$\lambda^{V}_{VBM}$ | -9.73 | -9.36 | -9.54 | 2.90 | 2.05 |
$\lambda^{V}_{CBM}$ | -7.82 | -7.67 | -7.76 | 4.68 | 3.93 |
|  |  |  |  |  |

Previous studies typically used the core state energy as the reference energy
level to correct the band edge energy. In pseudopotential method, the use of
pseudopotential to handle core electrons makes it impossible to obtain
accurate core state energy. Therefore, the deep valence (DV) state energy is
adopted as the reference energy level [36], and the DV state energy and the
band edge energy are obtained by the analyses of the density of states (DOS)
and the band structure along the high-symmetry points containing VBM and CBM,
respectively. Due to differences in calculation methods and parameters, the
DOS and band energies usually do not correspond perfectly. This increases the
uncertainty of the calculation results on the DP constants. In this work both
the DV state energy and the band edge energy are from the self-consistent
calculation with very dense _k_ points (30000 points in total). We take the
average energy of the lowest energy band as the DV state energy. Figs.
1(j)–1(l) show the DV state energy and the strain have a negative linear
relationship. The influence of the strain types on the DV state energy can be
almost ignored owing to small differences between the $\lambda_{dv}^{L}$,
$\lambda_{dv}^{S}$ and $\lambda_{dv}^{V}$ (-11.66 eV, -11.68 eV and -11.69
eV). More interestingly, the $\lambda_{dv}^{X}$ values under the compressive
and expanded strains are almost the same, and the difference between them is
less than 0.5 eV.

It is worth mentioning that we settle on the full potential rather than
pseudopotential method for electronic structure calculations and thus easily
obtain the DP constants of core states (1S and 2S) summarized in Table 2. The
core state energy like the DV state energy is insensitive to strain types (see
Fig.S1 of the Supporting Information), and the differences between them are
less than 0.02 eV.

Table 2: Summary of the calculated uniaxial/biaxial/volumetric DP constans (in units of eV) of core states (1S and 2S).  Core state | $\lambda_{core}^{L}$ | $\lambda_{core}^{S}$ | $\lambda_{core}^{V}$
---|---|---|---
1S | -12.4596 | -12.4592 | -12.4482
2S | -11.5886 | -11.5873 | -11.5954
|  |  |

After correcting the DP constants at the VBM and CBM using the core state
energy, we have found some abnormal conclusions as follows: (i) the DP
constants change from negative to positive values; (ii) the ratios of
$\lambda^{X}_{VBM}$ to $\lambda^{X}_{CBM}$ increase significantly; (iii) the
$|\lambda^{X}_{CBM}|$ values are much smaller than those of the
$|\lambda^{X}_{VBM}|$, which is the opposite situation before the correction,
because of the reference energy levels. In addition, different strain types
with the same strength lead to the difference that is comparable to, even
larger than the corrected DP constants. For instance, the difference between
the $\lambda^{S}_{VBM}$ and $\lambda^{L}_{VBM}$ is as large as 2.07 eV
significantly larger than the $\lambda^{L[1S]}_{VBM}$ and
$\lambda^{L[2S]}_{VBM}$. More importantly, choosing different core state
energies as the reference levels makes a difference to the DP constants,
leading to a large dissimilarity about 1 eV. Evidently, this means that the
strain has great influence on the deep core states.

Surprisedly, the corrected DP constants of the CBM, as well as the other
theoretical values, are in agreement well with the experimental data [28] that
was determined in the heavily As-doped Si with the n-type carrier density as
high as 5$\times 10^{21}cm^{-3}$. Unfortunately, the theoretical works neglect
an important experimental parameter of the carrier density related with the
Fermi energy level. Based on Boltzmann theory, the Fermi energy level
corresponds to the carrier density $n$ that is given by

$n=n_{0}-\int{f(\varepsilon,\varepsilon_{F},T)D(\varepsilon)d\varepsilon},$
(30)

where D, $\varepsilon_{F}$ and $\emph{n}_{0}$ are the total density of states,
the Fermi energy level and the number of the valence electrons in the unit
cell. Fig. 2 shows that the Fermi energy level of silicon at the fixed carrier
density is insensitive to temperature, however, is highly sensitive to the
carrier density at the fixed temperature. When the carrier density changes
from $2.7\times 10^{19}cm^{-3}$ to $5\times 10^{21}cm^{-3}$, the Fermi energy
level rises from 5.94 eV (corresponding to CBM) to 6.60 eV. This implies that
the experimental DP constants in Ref. [28] is for the Fermi energy level
rather than for the band edge where it is theoretically predicted in Refs.
[26, 27, 16]. Therefore, the high degree of agreement between theoretical and
experimental data strongly indicates that the theoretical predictions are
problematic. This is also confirmed by the fact that the calculated mobility
with the corrected DP constant severely deviates from the experimental data.

Figure 3: (Color online) (a) Plane-averaged potential along the (001)
orientation of silicon at V=$V_{0}$, 301 $V_{0}$. (b) The difference between
the plane-averaged potentials, (c) the difference between the plane-averaged
charge density at $V$ = $V_{0}$ and $V$ = 1.01 $V_{0}$ under volumetric
strain.

In this work we consider that the strain affects not only core electrons but
also valence electrons because the strain affects the entire periodic
potential especially near the silicon nucleus. Fig. 3(a) shows that both the
strength and shape of the periodic potential have undergone significant
changes in the entire region when the volume $V_{0}$ is expanded to 301
$V_{0}$. For instance, the periodic potential strength at the $V_{0}$ is lower
than -6.95 eV, however, it rises to roughly 0 eV at the 301 $V_{0}$, and the
potential shape change from the curve like the spring to an approximate
straight line. In fact, the two periodic potential values are not absolute due
to the different zero-potential points. However, the potential difference can
actually reflect the change in the periodic potential. The potential
difference in Fig. 3(b) shows the shape and strength of the periodic potential
have changed even if a small volume change of 1% occurs. Especially, the
closer to the silicon nucleus, the greater the periodic potential changes.

Table 3: The calculated effective mass of the silicon at different step sizes,
the corresponding mean values from averaging the last five step sizes, and the
experimental, theoretical data. In addition, the calculated and experimental
mobilities at 300 K are listed. The effective mass is in units of the electron
rest mass ($m_{0}$), the mobility is in units of $cm^{2}$/Vs, and the step
size is in units of 2$\pi$/4000$a_{0}$.

[b]  Step size $m_{l}$ $m_{t}$ $m_{hh}$ $m_{lh}$ 10 0.565 0.141 0.306 0.141 20
0.794 0.204 0.299 0.206 40 0.913 0.209 0.335 0.202 50 0.911 0.211 0.341 0.201
80 0.934 0.214 0.342 0.201 100 0.969 0.214 0.340 0.200 200 0.957 0.216 0.339
0.201 This work 0.937 0.213 0.339 0.201 Experiment 666The $m_{t}$ and $m_{l}$
values are taken from Ref. [37], the $m_{hh}$ and $m_{lh}$ values from [38].
0.9163 0.1905 0.43 0.19 Theory 2777Taken from Ref. [39]. 0.96 0.16 0.26 0.18
Temperature Calculated Measured $\mu_{p}$ $\mu_{n}$ $\mu_{p}$888Taken from
Ref. [40]. $\mu_{n}$999Taken from Ref. [41]. 300 K 706, 770 1771, 2631 360 to
510 1350 to 1750 99footnotetext: The $m_{t}$ and $m_{l}$ values are taken from
Ref. [37], the $m_{hh}$ and
$m_{lh}$ values from [38].

The electron density $\rho$(r) is determined by the probability density
related to the wave function that is independent of the zero-potential point.
The electron density difference truly reflects that the electron density
around the Si nucleus has more obvious change (see Fig. 3(c)). This change
originates from the atomic interaction change caused by volume variation, and
inevitably leads to the change for the core state energy. Generally, the
energy of the core electron is determined by the zero-potential point and the
atomic interactions composed of the exchange correction interaction and the
Coulomb interaction.

In this work we expand the cell volume large enough so that the effect of the
atomic interactions on the core electrons can be neglected. Obviously, the
change in the core state energy is mainly due to the zero-potential rather
than atomic interactions. Therefore, the correction for DP constant is
obtained easily.

Considering limited computing resources, we used the standard GGA rather than
mBJ-GGA strategy to perform the self-consistent calculation with rigorous
charge and energy convergence limits ($10^{-6}$ e and $10^{-6}$ eV). It is
worth noting to ensure the $R_{mt}\cdot K_{max}$ value in all calculations is
7.0 by setting suitable dimension parameters of WIEN2k code. Fig. 4 indicates
the relationship between the core state energy and the volumetric strains
seems to be linear. The slope values for 1S and 2S states are -0.024601 eV and
-0.024605 eV at 451 $V_{0}$, and decrease down to -0.021944 eV and -0.021987
eV at 501 $V_{0}$, respectively.

We consider the slight slope results from the zero-potential change induced by
the volumetric strain, and assume that the zero potential $E^{0}$ and the 1S
core state energy $E^{1S}$ are approximately equal, and have the property of
Coulomb potential

$E^{0}\propto\frac{k}{r}.$ (31)

Therefore, when the volume change from $V_{N}=NV_{0}$ to (1+$\delta^{V}$)
$V_{N}$, the zero potential difference $\Delta E^{0}_{N}$ can be described by

$\Delta
E^{0}_{N}\propto\frac{k}{\sqrt[3]{(1+\delta^{V})V_{N}}}-\frac{k}{\sqrt[3]{V_{N}}},$
(32)

So the correction for the DP constants at $V_{N}$ is about

$\lambda_{N}^{V}=\frac{\Delta E^{0}_{N}}{\delta^{V}}\approx\frac{\Delta
E^{1S}_{N}}{\delta^{V}}.$ (33)

And then one can obtain the correction for the DP constant at $V_{0}$
according the equation

$\lambda_{1}^{V}=\frac{\Delta
E^{0}_{1}}{\delta^{V}}\approx\frac{k}{\sqrt[3]{(1+\delta^{V})V_{0}}}-\frac{k}{\sqrt[3]{V_{0}}}\approx\sqrt[3]{N}\lambda^{V}_{N}$
(34)

Herein, taking N=501, the $\lambda_{1}^{V}$ is only about -0.1743 eV that is
insignificant compared to the $\lambda^{V}_{e(h)}$ without the correction. It
is worth mentioning that we also calculated the slope values at other large
volumes, and found that most of them follow this trend of the slope value
decreasing with the volume. However, a few slopes do not follow this trend,
implying that the uncertainty of zero-potential, but it is certain that the
corresponding $|\lambda_{1}^{V}|$ is not beyond 0.3 eV.

Figure 4: (Color online) Volumetric strain dependence of 1S and 2S core states
(a), (b) at $V$ = 501 $V_{0}$, and (c), (d) at $V$ = 451 $V_{0}$.

In addition, the effective mass has an important effect on the mobility, but
it is difficult to calculate because it is sensitive to the step size during
the numerical derivative process. The dependence of the effective mass with
step size for is summarised in Table 3. We average the last five values as the
effective mass values. The electron effective mass as well as the light-hole
mass is in good agreement well with measured and other calculated values. The
heavy-hole mass 0.34 $m_{0}$ is obviously larger than the other theoretical
values, however, it is much closer to the experimental value 0.43 $m_{0}$
[38]. We adopted the calculated elastic constant $c_{11}$ = 153 GPa that is
also in good agreement with experimental and theoretical values [42, 43, 44,
45]. Substituting the relevant parameters into the Eqs. (19) and (27) and
using the local DOS effective masses calculated by Eqs. (16) and (15)
respectively, the hole mobility at 300 K is calculated to be 704 cm2/Vs and
770 cm2/Vs that reasonably agrees with experimental value of 510 cm2/Vs [40]
(See Table 3). Adopting the mobility effective mass calculated by Eq. (29), we
obtain 1771 $cm^{2}$/Vs that agrees well with experimental data.

##  IV Conclusion

In summary, we improve the mobility formula based on the DP theory, and
combine with the first-principle method to calculate the electron and hole
mobilities for single crystal silicon. In the improved mobility formula, the
mobility effective mass is redefined, more importantly, a method is proposed
and applied to calculate the corrections for the DP constant. It is found that
the corrections for the DP are slight and thus can be neglected. Our data on
the mobility of the single crystal silicon is in good agreement with
experimental results, and thus the method can be effectively applied to the
prediction of the mobility in bulk materials.

## Acknowledgments

This work is supported by the Scientific and Technological project of Nanchang
Institute of Science and Technology (Grant No. NGKJ$-$22$-$09).

## References

  * Bardeen and Shockley [1950] J. Bardeen and W. Shockley, Phys. Rev. 80, 72 (1950).
  * Kartheuser and Rodriguez [1986] E. Kartheuser and S. Rodriguez, Physical Review B 33, 772 (1986).
  * Kaasbjerg, Thygesen, and Jacobsen [2012] K. Kaasbjerg, K. S. Thygesen, and K. W. Jacobsen, Physical Review B 85, 165440 (2012).
  * Murphy-Armando, Fagas, and Greer [2010] F. Murphy-Armando, G. Fagas, and J. Greer, Nano letters 10, 869 (2010).
  * Hong  _et al._ [2016] A. J. Hong, L. Li, R. He, J. J. Gong, Z. B. Yan, K. F. Wang, J. M. Liu, and Z. F. Ren, Scientific Reports 6, 22778 (2016).
  * Tang  _et al._ [2009] L. Tang, M. Long, D. Wang, and Z. Shuai, Science in China Series B: Chemistry 52, 1646 (2009).
  * Xi  _et al._ [2012] J. Xi, M. Long, L. Tang, D. Wang, and Z. Shuai, Nanoscale 4, 4348 (2012).
  * Kohn and Sham [1965] W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965).
  * Weinert, Wimmer, and Freeman [1982] M. Weinert, E. Wimmer, and A. J. Freeman, Phys. Rev. B 26, 4571 (1982).
  * Ihm, Zunger, and Cohen [1979] J. Ihm, A. Zunger, and M. L. Cohen, Journal of Physics C: Solid State Physics 12, 4409 (1979).
  * Janak [1974] J. F. Janak, Phys. Rev. B 9, 3985 (1974).
  * Janotti and Van de Walle [2007] A. Janotti and C. G. Van de Walle, Phys. Rev. B 75, 121201(R) (2007).
  * Van de Walle [1989] C. G. Van de Walle, Phys. Rev. B 39, 1871 (1989).
  * Van de Walle and Martin [1989] C. G. Van de Walle and R. M. Martin, Phys. Rev. Lett. 62, 2028 (1989).
  * Van de Walle and Neugebauer [1997] C. G. Van de Walle and J. Neugebauer, Applied Physics Letters 70, 2577 (1997).
  * Franceschetti, Wei, and Zunger [1994] A. Franceschetti, S.-H. Wei, and A. Zunger, Phys. Rev. B 50, 17797 (1994).
  * Wei and Zunger [1999] S.-H. Wei and A. Zunger, Phys. Rev. B 60, 5404 (1999).
  * Nolte, Walukiewicz, and Haller [1987] D. D. Nolte, W. Walukiewicz, and E. E. Haller, Phys. Rev. Lett. 59, 501 (1987).
  * Walukiewicz  _et al._ [1985] W. Walukiewicz, H. E. Ruda, J. Lagowski, and H. C. Gatos, Phys. Rev. B 32, 2645 (1985).
  * Lee  _et al._ [1979] H. J. Lee, J. Basinski, L. Y. Juravel, and J. C. Woolley, Canadian Journal of Physics 57, 233 (1979).
  * Pfeffer, Gorczyca, and Zawadzki [1984] P. Pfeffer, I. Gorczyca, and W. Zawadzki, Solid State Communications 51, 179 (1984).
  * Spitzer and Whelan [1959] W. G. Spitzer and J. M. Whelan, Phys. Rev. 114, 59 (1959).
  * Wang  _et al._ [2012] H. Wang, Y. Pei, A. D. LaLonde, and G. J. Snyder, Proceedings of the National Academy of Sciences 109, 9705 (2012).
  * Blacha, Presting, and Cardona [1984] A. Blacha, H. Presting, and M. Cardona, physica status solidi (b) 126, 11 (1984).
  * Vergés  _et al._ [1982] J. A. Vergés, D. Glötzel, M. Cardona, and O. K. Andersen, physica status solidi (b) 113, 519 (1982).
  * Li, Gong, and Wei [2006] Y.-H. Li, X. G. Gong, and S.-H. Wei, Applied Physics Letters 88 (2006), 10.1063/1.2168254, 042104.
  * Resta, Colombo, and Baroni [1990] R. Resta, L. Colombo, and S. Baroni, Phys. Rev. B 41, 12358 (1990).
  * Cargill, Angilello, and Kavanagh [1988] G. S. Cargill, J. Angilello, and K. L. Kavanagh, Phys. Rev. Lett. 61, 1748 (1988).
  * Luo  _et al._ [2014] Y. Luo, J. Wang, J. Wang, J. Li, and Z. Hu, Journal of the American Ceramic Society 97, 945 (2014).
  * [30] X. Liu, T. Zhu, H. Wang, L. Hu, H. Xie, G. Jiang, G. J. Snyder, and X. Zhao, Advanced Energy Materials 3, 1238, https://onlinelibrary.wiley.com/doi/pdf/10.1002/aenm.201300174 .
  * Batchelder and Simmons [2004] D. N. Batchelder and R. O. Simmons, The Journal of Chemical Physics 41, 2324 (2004).
  * Bludau, Onton, and Heinke [2003] W. Bludau, A. Onton, and W. Heinke, Journal of Applied Physics 45, 1846 (2003).
  * Tran and Blaha [2009] F. Tran and P. Blaha, Phys. Rev. Lett. 102, 226401 (2009).
  * Schwarz, Blaha, and Madsen [2002] K. Schwarz, P. Blaha, and G. Madsen, Computer Physics Communications 147, 71 (2002).
  * Camargo-Martínez and Baquero [2012] J. A. Camargo-Martínez and R. Baquero, Phys. Rev. B 86, 195106 (2012).
  * Zhao  _et al._ [2016] T. Zhao, W. Shi, J. Xi, D. Wang, and Z. Shuai, Scientific Reports 6, 19968 (2016).
  * Hensel, Hasegawa, and Nakayama [1965] J. C. Hensel, H. Hasegawa, and M. Nakayama, Phys. Rev. 138, A225 (1965).
  * Dexter, Zeiger, and Lax [1956] R. N. Dexter, H. J. Zeiger, and B. Lax, Phys. Rev. 104, 637 (1956).
  * Ramos  _et al._ [2001] L. E. Ramos, L. K. Teles, L. M. R. Scolfaro, J. L. P. Castineira, A. L. Rosa, and J. R. Leite, Phys. Rev. B 63, 165210 (2001).
  * Cronemeyer [1957] D. C. Cronemeyer, Phys. Rev. 105, 522 (1957).
  * Rode [1972] D. L. Rode, physica status solidi (b) 53, 245 (1972).
  * Łopuszyński and Majewski [2007] M. Łopuszyński and J. A. Majewski, Phys. Rev. B 76, 045202 (2007).
  * McSkimin, H. J. and Andreatch, P., Jr. [2004] McSkimin, H. J. and Andreatch, P., Jr., Journal of Applied Physics 35, 3312 (2004).
  * Zhao, Winey, and Gupta [2007] J. Zhao, J. M. Winey, and Y. M. Gupta, Phys. Rev. B 75, 094105 (2007).
  * Nielsen and Martin [1985] O. H. Nielsen and R. M. Martin, Phys. Rev. B 32, 3792 (1985).