# Joint Blocklength and Location Optimization for URLLC-enabled UAV Relay
Systems

Cunhua Pan, Hong Ren, Yansha Deng, Maged Elkashlan, and Arumugam Nallanathan

###### Abstract

This letter considers the unmanned aerial vehicle (UAV)-enabled relay system
to deliver command information under ultra-reliable and low-latency
communication (URLLC) requirements. We aim to jointly optimize the blocklength
allocation and the UAV’s location to minimize the decoding error probability
subject to the latency requirement. The achievable data rate under finite
blocklength regime is adopted. A novel perturbation-based iterative algorithm
is proposed to solve this problem. Simulation results show that the proposed
algorithm can achieve the same performance as the exhaustive search method,
and significantly outperforms the existing algorithms.

##  I Introduction

UAV-assisted communication has attracted extensive attention due to its fast
deployment and favorable channel gain [1]. UAVs can also serve as relays to
provide wireless connectivity between two devices without direct communication
links [1, 2, 3, 4, 5]. Joint relay trajectory and power allocation was studied
in [2]. In [3], UAV node placement and communication resource allocation were
jointly optimized. In [4], Zhang _et al._ studied the joint trajectory and
power optimization to minimize the outage probability. In [5], the throughput
maximization problem was studied for a two-user broadcast channel, which can
be regarded as a decode-forward relay system where each hop has the same rate.

In 2017, URLLC has been regarded as one of three pillar applications that
should be supported in the 5G communications [6]. Applications requiring URLLC
services include factory automation, autonomous driving, remote surgery, etc.
In URLLC, short packet transmission is normally selected to support the low-
latency transmission [7]. In this case, the conventional Shannon’s capacity
based on the law of large numbers is no longer applicable. The achievable
capacity under short packet regime was first derived in [8], which is a
complicated function of the system parameters.

Recently, [9] and [10] considered the delay issues in UAV communications. Mean
packet transmission delay minimization problem was studied in [9] with multi-
layer UAVs, while minimum-rate ratio for each user was considered in [10] to
flexibly adjust the percentage of its delay-constrained data traffic. However,
the latency requirement was not considered in [9, 10], and the Shannon’s
capacity formula was adopted.

In this paper, we consider a downlink communication system in a frontline as
shown in Fig. 1, where a central controller needs to send command information
to a distant robot that performs certain reconnaissance missions in a military
area. For concealment, shelters with thick cement/metal walls are built
between the outside controller and the military area. Hence, the channel gain
between the controller and the robot is weak and negligible, and requires a
UAV to fly above the shelter to assist the transmission between the controller
and the robot. We aim for jointly optimizing the location of the UAV and
blocklength allocation of the UAV and the controller to minimize the decoding
error probability subject to the latency and location constraints. To solve
this problem, we propose a novel perturbation-based iterative algorithm to
alternatively optimize the location and the blocklength. Our results show that
the proposed low-complexity algorithm achieves almost the same performance as
the exhaustive search method, and performs much better than the existing
algorithms.

##  II System Model

As shown in Fig. 1, we consider a two-dimensional UAV-enabled military
surveillance scenario111In this paper, we consider only one UAV and one robot
since it may be easily discovered by the enemies with more UAVs and robots. In
addition, the expression of the decoding error probability for multiple UAVs
and robots is very complicated and difficult to optimize, which will be left
for future work., where the UAV hovers at $(x,H)$ above the horizontal line
between the controller and the robot222It is obvious that better system
performance can be achieved with high channel power gains, which is the case
when the UAV hovers above the line between the controller and the robot., with
$H$ as the fixed altitude. The locations of controller and robot are $(0,0)$
and $(D,0)$. The packet size of the command signal is $L$ bits, whose
transmission needs to be completed within $T_{\rm{max}}$ seconds. Then, the
overall blocklength is $M=BT_{\rm{max}}$ [7], where $B$ is the system
bandwidth. Each transmission period has two phases, i.e., the first phase
corresponds to the transmission from controller to UAV, while the second is
from UAV to robot. The blocklength allocated for each phase is given by
$m_{1}$ and $m_{2}$, respectively. The transmission powers from controller and
UAV are fixed as $P_{1}$ and $P_{2}$, respectively.

Figure 1: UAV relay system for delivering URLLC services.

The channel power gain from the controller to the UAV, and that from the UAV
to the robot are denoted as ${h_{1}}$ and ${h_{2}}$, respectively. According
to the measurement result in [11], the LOS probability is close to one when
the UAV is above a certain altitude (e.g. 120 m), and the free space channel
model can be adopted. Thus, $h_{1}$ and $h_{2}$ can be represented as

${h_{1}}=\frac{{{\beta_{0}}}}{{{H^{2}}+{x^{2}}}},{h_{2}}=\frac{{{\beta_{0}}}}{{{H^{2}}+{{(D-x)}^{2}}}},\vspace{-0.1cm}$
(1)

where $\beta_{0}$ is channel power gain at a reference distance of $d_{0}=1$
meter.

According to [8], to transmit a short packet of size $L$ within $m_{1}$
symbols, the decoding error at the UAV is given by
${\varepsilon_{1}}=Q\left({f\left({{\gamma_{1}},{m_{1}},L}\right)}\right)$,
where $f\left({\gamma_{1},m_{1},L}\right)=\ln
2\sqrt{\frac{m_{1}}{V_{1}}}\left({{{\log}_{2}}(1+\gamma_{1})-\frac{L}{m_{1}}}\right)$
with $\gamma_{1}=P_{1}h_{1}$ 333The noise is normalized to unit. and
$V_{1}=1-(1+\gamma_{1})^{-2}$. Similarly, the decoding error probability at
the robot is given by
${\varepsilon_{2}}=Q\left({f\left({{\gamma_{2}},{m_{2}},L}\right)}\right)$,
where $\gamma_{2}=P_{2}h_{2}$.

We consider that the UAV acts as a decode-and-forward (DF) relay. Then, the
overall decoding error probability from the controller to the robot is given
by

$\varepsilon={\varepsilon_{1}}+(1-{\varepsilon_{1}}){\varepsilon_{2}}.\vspace{-0.2cm}$
(2)

To enable URLLC, we aim to jointly optimize the location of the UAV and the
blocklength for two phases to minimize the overall decoding error probability
under the latency/blocklength constraint. Thus, the optimization problem can
be formulated as

$\displaystyle\mathop{\min}\limits_{\left\\{{{m_{1}},{m_{2}},x}\right\\}}\;\;\;$
$\displaystyle{{\varepsilon}}$ (3a) $\displaystyle{\rm{s.t.}}\;\;\;$
$\displaystyle d_{1}\leq x\leq d_{2},$ (3b) $\displaystyle m_{1}+m_{2}=M,$
(3c) $\displaystyle m_{1},m_{2}\in\mathbb{Z},\vspace{-0.1cm}$ (3d)

where (3b) specifies the feasible region for $x$ that mainly depends on the
shape of the shelters, $\mathbb{Z}$ is the positive integer set. Problem (3)
is difficult to solve due to: 1) the objective function is not jointly convex
w.r.t. the optimization variables; 2) even with fixed blocklength or location,
the objective function is not convex w.r.t. the location or the blocklength.

##  III Low-complexity Algorithm

In this section, we develop a low-complexity iterative algorithm to solve
Problem (3). In particular, we first find the optimized blocklength with fixed
location $x$, and then find the optimal location by fixing the blocklength
allocation, and at last alternatively solve each subproblem until convergence.

###  III-A Optimize Blocklength Allocation with Fixed $x$

With fixed UAV location $x$, the channel gains $h_{1}$ and $h_{2}$ are fixed
according to (1). Then, we only need to optimize the blocklength. Remind that
the expression of objective function in (2) is very complicated. To handle
this difficulty, we note that to guarantee the extremely low error probability
for the whole link, the error probability for each link should be sufficiently
small. In this case, we can approximate the overall error probability as
$\varepsilon(m_{1})\approx{\varepsilon_{1}}(m_{1})+{\varepsilon_{2}}(m_{1})\triangleq\tilde{\varepsilon}(m_{1})$.
In the following theorem, we prove that $\tilde{\varepsilon}(m_{1})$ is a
convex function of $m_{1}$.

_Theorem 1_ : We assume that $m_{1}$ is a continuous variable. Given the
channel gains $h_{1}$ and $h_{2}$, $\tilde{\varepsilon}(m_{1})$ is a convex
function w.r.t. $m_{1}$.

_Proof_ : We first prove that ${\varepsilon_{1}}(m_{1})$ is a convex function
w.r.t. $m_{1}$, then the convexity of function ${\varepsilon_{2}}(m_{1})$ can
be proved in a similar method by replacing $m_{2}$ with $M-m_{1}$. Thus, the
summation of two convex functions $\tilde{\varepsilon}(m_{1})$ is a convex
function.

To simplify the notation, we denote
$f_{1}(m_{1})={f\left({{\gamma_{1}},{m_{1}},L}\right)}$. The first order and
second order of ${\varepsilon_{1}}(m_{1})$ w.r.t. $m_{1}$ are given by

$\displaystyle\varepsilon_{1}^{\prime}(m_{1})\\!\\!\\!\\!\\!\\!$
$\displaystyle=$
$\displaystyle\\!\\!\\!\\!\\!-\frac{1}{{\sqrt{2\pi}}}{e^{-\frac{{f_{1}^{2}(m_{1})}}{2}}}{f_{1}^{\prime}(m_{1})},$
(4) $\displaystyle\varepsilon_{1}^{\prime\prime}(m_{1})\\!\\!\\!\\!\\!\\!$
$\displaystyle=$
$\displaystyle\\!\\!\\!\\!\\!\\!\frac{1}{{\sqrt{2\pi}}}{e^{-\frac{{f_{1}^{2}(m_{1})}}{2}}}\\!\left({{f_{1}(m_{1})}{{\left({{f_{1}^{\prime}(m_{1})}}\right)}^{2}}\\!\\!-\\!\\!f_{1}^{\prime\prime}(m_{1})}\right).$
(5)

Since the first term in the brackets of (5) is positive, we only need to prove
${f_{1}^{\prime\prime}}(m_{1})\leq 0$ so that ${\varepsilon_{1}}(m_{1})$ is a
convex function.

To this end, we first simplify the function $f_{1}$ as

${f_{1}}(m_{1})={A_{1}}\sqrt{{m_{1}}}\left({{C_{1}}-\frac{L}{{{m_{1}}}}}\right),$
(6)

where $A_{1}$ and $C_{1}$ are constants given by ${A_{1}}={{\ln
2}\mathord{\left/{\vphantom{{\ln
2}{\sqrt{{V_{1}}}}}}\right.\kern-1.2pt}{\sqrt{{V_{1}}}}}$ and
${C_{1}}={\log_{2}}(1+{\gamma_{1}})$, respectively.

The first and second derivative of ${f_{1}}(m_{1})$ w.r.t. $m_{1}$ are

$\displaystyle{f_{1}^{\prime}}(m_{1})$ $\displaystyle=$
$\displaystyle\frac{1}{2}{A_{1}}{C_{1}}m_{1}^{-\frac{1}{2}}+\frac{1}{2}{A_{1}}Lm_{1}^{-\frac{3}{2}},$
(7) $\displaystyle{f_{1}^{\prime\prime}}(m_{1})$ $\displaystyle=$
$\displaystyle-\frac{1}{4}{A_{1}}{C_{1}}m_{1}^{-\frac{3}{2}}-\frac{3}{4}{A_{1}}Lm_{1}^{-\frac{5}{2}}<0,$
(8)

which show that ${f_{1}}(m_{1})$ is a strictly concave function of $m_{1}$.
Then ${\varepsilon_{1}}(m_{1})$ is a convex function w.r.t. $m_{1}$. By using
the similar derivations, we show that ${\varepsilon_{2}}(m_{1})$ is also a
convex function w.r.t. $m_{1}$. This completes the proof.

Based on Theorem 1, we apply the bisection search method to find the optimal
blocklength allocation through solving the following equation:

$\tilde{\varepsilon}^{\prime}(m_{1})=\varepsilon_{1}^{\prime}(m_{1})+\varepsilon_{2}^{\prime}(m_{1})=0.$
(9)

The details of searching for the optimal blocklength allocation are provided
in Algorithm 1.

1 Initialize $m_{1}^{\rm{lb}}=1,m_{1}^{\rm{ub}}=M-1$ and the error tolerance
$\delta=0.5$;

2while  _ $m_{1}^{\rm{ub}}-m_{1}^{\rm{lb}}>\delta$ _ do

3 Set $m_{1}^{\rm{mid}}=(m_{1}^{\rm{lb}}+m_{1}^{\rm{ub}})/2$.

4 if
_${\left.{\tilde{\varepsilon}^{\prime}(m_{1})}\right|_{{m_{1}}=m_{1}^{{\rm{mid}}}}}
>0$_ then

5 Set $m_{1}^{\rm{ub}}=m_{1}^{{\rm{mid}}}$.

6 else

7 Set $m_{1}^{\rm{lb}}=m_{1}^{{\rm{mid}}}$.

8 end if

9

10 end while

Return the optimal $m_{1}$ as
$m_{1}^{*}=\mathop{\arg\min}\limits_{\left\\{{\left\lfloor{m_{1}^{{\rm{mid}}}}\right\rfloor,\left\lceil{m_{1}^{{\rm{mid}}}}\right\rceil}\right\\}}\tilde{\varepsilon}$
and the optimal $m_{2}$ as $m_{2}^{*}=M-m_{1}^{*}$.

Algorithm 1 Find the Optimal $m_{1}$ and $m_{2}$ with Fixed $x$

###  III-B Optimal Location Optimization with Fixed $m_{1}$ and $m_{2}$

In this subsection, we aim for optimizing location $x$ with given $m_{1}$ and
$m_{2}$. However, the objective function $\tilde{\varepsilon}(x)$ is not a
convex function w.r.t. $x$. In the following, we numerically show that
$\tilde{\varepsilon}(x)$ has only one local minimum point. Hence, there only
exists only one solution that minimizes $\tilde{\varepsilon}(x)$.

To provide clear explanations, we define a new function
$g(x)\buildrel\Delta\over{=}{\ln}(\tilde{\varepsilon}(x))$, which has the same
monotonic property of $\tilde{\varepsilon}(x)$. We first obtain the first and
second derivative of $g(x)$ w.r.t. $x$ as follows. We define
$f_{1}(x)={f\left({{\gamma_{1}}(x),{m_{1}},L}\right)}$ and
$f_{2}(x)={f\left({{\gamma_{2}}(x),{m_{2}},L}\right)}$, where
${\gamma_{1}}(x)$ and ${\gamma_{2}}(x)$ are given by
${\gamma_{1}}(x)=\frac{{{P_{1}\beta_{0}}}}{{{H^{2}}+{x^{2}}}}$ and
${\gamma_{2}}(x)=\frac{{{P_{2}\beta_{0}}}}{{{H^{2}}+{(D-x)^{2}}}}$,
respectively.

The first derivative of $g(x)$ w.r.t. $x$ is given by

$g^{\prime}(x)=\frac{{{{\varepsilon^{\prime}_{1}}}(x)+{{\varepsilon^{\prime}_{2}}}(x)}}{{{\varepsilon_{1}(x)}+{\varepsilon_{2}(x)}}},$
(10)

where ${{\varepsilon^{\prime}_{i}}}(x)$ is given by

${{\varepsilon^{\prime}_{i}}}(x)=\frac{{\partial{\varepsilon_{i}}({\gamma_{i}})}}{{\partial{\gamma_{i}}}}{{\gamma^{\prime}_{i}}}(x),i=1,2,$
(11)

with ${{\gamma^{\prime}_{i}}}(x)$ given by

${{\gamma^{\prime}_{1}}}(x)=-\frac{{2{P_{1}}{\beta_{0}}x}}{{{{\left({{H^{2}}+{x^{2}}}\right)}^{2}}}},{{\gamma^{\prime}_{2}}}(x)=\frac{{2{P_{2}}{\beta_{0}}(D-x)}}{{{{\left({{H^{2}}+{{(D-x)}^{2}}}\right)}^{2}}}},$
(12)

and $\frac{{\partial{\varepsilon_{i}}({\gamma_{i}})}}{{\partial{\gamma_{i}}}}$
given by

$\frac{{\partial{\varepsilon_{i}}({\gamma_{i}})}}{{\partial{\gamma_{i}}}}=-\frac{1}{{\sqrt{2\pi}}}{e^{-\frac{{f_{i}^{2}({\gamma_{i}})}}{2}}}\frac{{\partial{f_{i}}({\gamma_{i}})}}{{\partial{\gamma_{i}}}}.$
(13)

In (13), $\frac{{\partial{f_{i}}({\gamma_{i}})}}{{\partial{\gamma_{i}}}}$ is

$\vspace{-0.3cm}\frac{{\partial{f_{i}}({\gamma_{i}})}}{{\partial{\gamma_{i}}}}=\sqrt{{m_{i}}}\frac{{1-\ln
2\frac{{{{\log}_{2}}(1+{\gamma_{i}})-\frac{L}{{{m_{i}}}}}}{{{{(1+{\gamma_{i}})}^{2}}-1}}}}{{\sqrt{{{(1+{\gamma_{i}})}^{2}}-1}}}.$
(14)

The second derivative of $g(x)$ w.r.t. $x$ can be calculated as

$g^{\prime\prime}(x)\\!=\\!\frac{{\left({{{\varepsilon^{\prime\prime}_{1}}}(x)+{{\varepsilon^{\prime\prime}_{2}}}(x)}\right)\left({{\varepsilon_{1}}(x)+{\varepsilon_{2}}(x)}\right)-{{\left({{{\varepsilon^{\prime}_{1}}}(x)+{{\varepsilon^{\prime}_{2}}}(x)}\\!\right)}^{2}}}}{{{{\left({{\varepsilon_{1}}(x)+{\varepsilon_{2}}(x)}\right)}^{2}}}},$

where ${{{\varepsilon^{\prime\prime}_{i}}}(x)},i=1,2$ are given by

${{\varepsilon^{\prime\prime}_{i}}}(x)=\frac{{{\partial^{2}}{\varepsilon_{i}}({\gamma_{i}})}}{{\partial\gamma_{i}^{2}}}{\left({{{\gamma^{\prime}_{i}}}(x)}\right)^{2}}+\frac{{\partial{\varepsilon_{i}}({\gamma_{i}})}}{{\partial{\gamma_{i}}}}{{\gamma^{\prime\prime}_{i}}}(x),$
(15)

with ${{\gamma^{\prime\prime}_{i}}}(x)$ given by

$\displaystyle{{\gamma^{\prime\prime}_{1}}}(x)\\!\\!\\!\\!\\!$
$\displaystyle=$
$\displaystyle\\!\\!\\!\\!\\!\frac{{6{P_{1}}{\beta_{0}}{x^{4}}+4{P_{1}}{\beta_{0}}{H^{2}}{x^{2}}-2{P_{1}}{\beta_{0}}{H^{4}}}}{{{{\left({{H^{2}}+{x^{2}}}\right)}^{4}}}},$
$\displaystyle{{\gamma^{\prime\prime}_{2}}}(x)\\!\\!\\!\\!\\!$
$\displaystyle=$
$\displaystyle\\!\\!\\!\\!\\!\frac{{6{P_{2}}{\beta_{0}}{{(D\\!-\\!x)}^{4}}\\!\\!+\\!4{P_{2}}{\beta_{0}}{H^{2}}{{(D\\!-\\!x)}^{2}}\\!\\!-\\!\\!2{P_{2}}{\beta_{0}}{H^{4}}}}{{{{\left(\\!{{H^{2}}+{{(D\\!-\\!x)}^{2}}}\\!\right)}^{4}}}},$

and
$\frac{{{\partial^{2}}{\varepsilon_{i}}({\gamma_{i}})}}{{\partial\gamma_{i}^{2}}}$
given by

$\\!\\!\frac{{{\partial^{2}}{\varepsilon_{i}}({\gamma_{i}})}}{{\partial\gamma_{i}^{2}}}\\!\\!=\\!\\!\frac{1}{{\sqrt{2\pi}}}{e^{-\frac{{f_{i}^{2}({\gamma_{i}})}}{2}}}\\!\\!\left(\\!\\!{{f_{i}}({\gamma_{i}})\\!{{\left(\\!{\frac{{\partial{f_{i}}({\gamma_{i}})}}{{\partial{\gamma_{i}}}}}\\!\right)}^{2}}\\!\\!\\!-\\!\frac{{{\partial^{2}}{f_{i}}({\gamma_{i}})}}{{\partial\gamma_{i}^{2}}}}\\!\\!\right).$
(16)

In (16),
${\frac{{{\partial^{2}}{f_{i}}({\gamma_{i}})}}{{\partial\gamma_{i}^{2}}}}$ is
given by

$\displaystyle\frac{{{\partial^{2}}{f_{i}}({\gamma_{i}})}}{{\partial\gamma_{i}^{2}}}\\!\\!$
$\displaystyle=$
$\displaystyle\\!\\!\\!\\!c_{i}\left({-\frac{1}{{1+{\gamma_{i}}}}-(1+{\gamma_{i}})}\right)\left({{{\left({1+{\gamma_{i}}}\right)}^{2}}-1}\right)\\!+$
(17) $\displaystyle
3c_{i}(1+{\gamma_{i}})\left({{{\log}_{2}}(1+{\gamma_{i}})\\!-\\!\frac{L}{{{m_{i}}}}}\right)\ln
2,$

where
$c_{i}={{\sqrt{{m_{i}}}}\mathord{\left/{\vphantom{{\sqrt{{m_{i}}}}{{{\left({{{\left({1+{\gamma_{i}}}\right)}^{2}}-1}\right)}^{\frac{5}{2}}}}}}\right.\kern-1.2pt}{{{\left({{{\left({1+{\gamma_{i}}}\right)}^{2}}-1}\right)}^{5/2}}}}$.

Figure 2: $g(x)$, $g^{\prime}(x)$ and $g^{\prime\prime}(x)$ versus $x$, where
the system parameters are the same as those in the simulation section.

Unfortunately, $g(x)$ is not a convex function w.r.t. $x$ in the whole region
of $x$. Fig. 2 plots the functions $g(x)$, $g^{\prime}(x)$ and
$g^{\prime\prime}(x)$ versus $x$. It is observed from this figure that when
$0\rm{m}<x<142.5\rm{m}$, $g^{\prime\prime}(x)<0$, which means $g(x)$ is a
concave function. On the other hand, when $142.5\rm{m}<x<186.5\rm{m}$,
$g^{\prime\prime}(x)>0$ and $g(x)$ is a convex function. Finally, when
$186.5\rm{m}<x<200\rm{m}$, $g(x)$ becomes a concave function again. However,
it is observed that $g^{\prime}(x)>0$ when $0\rm{m}<x<163.4\rm{m}$, and
$g^{\prime}(x)<0$ when $163.4\rm{m}<x<200\rm{m}$. This means function $g(x)$
first decreases with $x$ for $0\rm{m}<x<163.4\rm{m}$ and then increases with
$x$ for $163.4\rm{m}<x<200\rm{m}$, and there exists only one minimum value. It
is difficult to rigorously prove this due to the cumbersome expression of
$g^{\prime}(x)$. We check all the other simulation parameters, all numerical
results show the same trend of function $g(x)$ 444The rational behind this is
that when x is small, the link from UAV to the robot will be the bottleneck of
the whole link, and vice versa. . As a result, the bisection search can be
used to find the root of $g^{\prime}(x)$, details of which can be found in the
following algorithm.

1 Initialize $x^{\rm{lb}}=d_{1},x^{\rm{ub}}=d_{2}$ and error tolerance
$\zeta=0.1$;

2if  _${\left.{g^{\prime}(x)}\right|_{x={d_{2}}}} <0$_ then

3 Return the optimal $x^{*}=d_{2}$.

4else

5 if  _${\left.{g^{\prime}(x)}\right|_{x={d_{1}}}} >0$_ then

6 Return the optimal $x^{*}=d_{1}$.

7 else

8 while  _$x^{\rm{ub}}-x^{\rm{lb}} >\zeta$_ do

9 Set $x^{\rm{mid}}=(x^{\rm{ub}}+x^{\rm{lb}})/2$.

10 if  _${\left.{g^{\prime}(x)}\right|_{x=x^{\rm{mid}}}} >0$_ then

11 Set $x^{\rm{up}}=x^{\rm{mid}}$.

12 else

13 Set $x^{\rm{lb}}=x^{\rm{mid}}$.

14 end if

15

16 end while

17

18 end if

19

20 end if

Algorithm 2 Find the Optimal $x$ with Fixed $m_{1}$ and $m_{2}$

###  III-C Overall Algorithm and Analysis

1 Initialize $m_{1}^{(0)},m_{2}^{(0)}$ and $x^{(0)}$, integer parameter
$N_{\rm{max}}$, iterative index $t=1$, maximum iterative times $t_{\rm{max}}$;

2repeat

3 With given $x^{(t-1)}$, obtain the optimal $m_{1}^{(t)}$ by using Algorithm
1;

4 Generate two random integer values $n^{l}$ and $n^{r}$ within
$[1,N_{\rm{max}}]$. Use Algorithm 2 to calculate the optimal location when
$m_{1}=m_{1}^{(t)}-n^{l},m_{1}^{(t)},$ $m_{1}^{(t)}+n^{r}$. The corresponding
optimal $x$ is denoted as ${\cal X}^{(t)}=\\{x_{l}^{(t)},$
$x_{m}^{(t)},x_{r}^{(t)}\\}$. Choose the optimal $x$ from ${\cal X}^{(t)}$
with the minimum $\tilde{\varepsilon}$, and denote it as $x^{(t)}$.

5 Set $t=t+1$;

6until  _ $t\geq t_{\rm{max}}$ _;

Algorithm 3 Iterative Algorithm for Solving Problem (3)

It is noted from simulations that the conventional block coordinate descent
method, which directly iterates between blocklength and location, is very
likely to get stuck at the initial point. To overcome this issue, we introduce
a small perturbation for the blocklength $m_{1}$ in each iteration as shown in
line 4 of Algorithm 3, where $N_{\rm{max}}$ is a small integer.

The overall algorithm for solving Problem (3) is provided in Algorithm 3. The
convergence of this algorithm is guaranteed since the objective value
decreases in each step and the value is lower-bounded by zero. The complexity
of this algorithm is analyzed as follows. In each iteration of Algorithm 3, we
need to run both Algorithm 1 only once and Algorithm 2 for three times. The
complexity of Algorithm 1 is given by $Q_{1}={\cal
O}\left({{{\log}_{2}}(M/\delta)}\right)$, while that of Algorithm 2 is
$Q_{2}={\cal
O}\left({{{\log}_{2}}\left({{{({d_{2}}-{d_{1}})}\mathord{\left/{\vphantom{{({d_{2}}-{d_{1}})}\zeta}}\right.\kern-1.2pt}\zeta}}\right)}\right)$.
Hence, the overall complexity of Algorithm 3 is given by ${\cal
O}\left({{n_{\max}}({Q_{1}}+{Q_{2}})}\right)$. In simulations, the algorithm
generally converges within ten iterations.

By using the similar analysis, the complexity of exhaustive search method is
given by ${\cal O}\left(M(d_{2}-d_{1})/\zeta\right)$, which is significantly
higher that of Algorithm 3.

Since the original problem in (3) is non-convex, Algorithm 3 cannot be
guaranteed to yield a globally optimal solution. However, from the simulation
results, we can find that this algorithm can achieve the same performance as
the exhaustive search method.

##  IV Simulation Results

We now perform simulation results to show the performance of our proposed
algorithm. The system parameters are set as follows: system bandwidth of $B=1\
{\rm{MHz}}$, $D=200$ m, $H=120$ m, $d_{1}=30$ m, $d_{2}=130$ m, $L=100$ bits,
$M=100$, $\beta_{0}=50$ dB, $P_{1}=3$ Watt, $P_{2}=1$ Watt, $N_{\rm{max}}=3$,
$t_{\rm{max}}=10$. The system transmission delay duration is set as
$T_{\rm{max}}=100\ {\rm{us}}$. Thus, the total number of symbols is
$M=BT_{\rm{max}}=100$.

Figure 3: Convergence behaviour of Algorithm 3 for various $H$.

In Fig. 3, we plot the convergence behaviour of Algorithm 3 for various $H$.
It is shown in Fig. 3 that the algorithm converges rapidly and generally ten
iterations are enough for convergence for all considered $H$, which indicates
that our algorithm has a low complexity.

Figure 4: Performance comparison for various algorithms.

In Fig. 4, we compared the proposed algorithm with the following algorithms:
1) exhaustive search algorithm (labeled as ‘Exhaus.’), 2) optimal blocklength
allocation with fixed location $x=(d_{1}+d_{2})/2$ (labeled as ‘Fixedx’), and
3) the optimal location with fixed blocklength $m_{1}=M/2$ (labeled as
‘Fixedm’). It is seen from Fig. 4 that the proposed algorithm achieves the
same performance as that of the exhaustive search algorithm, and significantly
outperforms the other two algorithms, which emphasizes the importance of joint
optimization.

##  V Conclusions

This paper studied the joint location and blocklength allocation for UAV relay
system with URLLC requirement. An effective low-complexity iterative algorithm
was proposed to solve the optimization problem. Each subproblem can be solved
by using the bisection search method. Simulation results confirm that the
proposed algorithm achieved the same performance as that of the exhaustive
search method, and has superior performance over the existing algorithms.

## References

  * [1] Y. Zeng, R. Zhang, and T. J. Lim, “Wireless communications with unmanned aerial vehicles: opportunities and challenges,” _IEEE Commun. Mag.,_ , vol. 54, no. 5, pp. 36–42, May 2016.
  * [2] ——, “Throughput maximization for UAV-enabled mobile relaying systems,” _IEEE Trans. Commun._ , vol. 64, no. 12, pp. 4983–4996, Dec 2016.
  * [3] R. Fan, J. Cui, S. Jin _et al._ , “Optimal node placement and resource allocation for UAV relaying network,” _IEEE Commun. Lett._ , vol. 22, no. 4, pp. 808–811, April 2018.
  * [4] S. Zhang, H. Zhang _et al._ , “Joint trajectory and power optimization for UAV relay networks,” _IEEE Commun. Lett._ , vol. 22, no. 1, pp. 161–164, Jan 2018.
  * [5] Q. Wu, J. Xu, and R. Zhang, “Capacity characterization of UAV-enabled two-user broadcast channel,” _IEEE J. Sel. Areas Commun._ , vol. 36, no. 9, pp. 1955–1971, Sep. 2018.
  * [6] M. Shafi, A. F. Molisch _et al._ , “5G: A tutorial overview of standards, trials, challenges, deployment, and practice,” _IEEE J. Sel. Areas Commun._ , vol. 35, no. 6, pp. 1201–1221, June 2017.
  * [7] G. Durisi, T. Koch, and P. Popovski, “Toward massive, ultrareliable, and low-latency wireless communication with short packets,” _Proc. IEEE_ , vol. 104, no. 9, pp. 1711–1726, Sept 2016.
  * [8] Y. Polyanskiy, H. V. Poor, and S. Verdu, “Channel coding rate in the finite blocklength regime,” _IEEE Trans. Inf. Theory_ , vol. 56, no. 5, pp. 2307–2359, May 2010.
  * [9] J. Li and Y. Han, “Optimal resource allocation for packet delay minimization in multi-layer UAV networks,” _IEEE Commun. Lett._ , vol. 21, no. 3, pp. 580–583, March 2017.
  * [10] Q. Wu and R. Zhang, “Common throughput maximization in UAV-enabled OFDMA systems with delay consideration,” _IEEE Trans. Commun._ , vol. 66, no. 12, pp. 6614–6627, Dec 2018.
  * [11] X. Lin, V. Yajnanarayana, S. D. Muruganathan, S. Gao, H. Asplund, H. Maattanen, M. Bergstrom, S. Euler, and Y. . E. Wang, “The sky is not the limit: Lte for unmanned aerial vehicles,” _IEEE Commun. Mag._ , vol. 56, no. 4, pp. 204–210, April 2018.


††thanks: Corresponding author
<EMAIL_ADDRESS>Corresponding author
<EMAIL_ADDRESS>
# Radically Uniform Spike Trains in Optically Injected Quantum Cascade
Oscillators

Yibo Peng  School of Information Science and Technology, ShanghaiTech
University, Shanghai 201210, China  Shanghai Institute of Microsystem and
Information Technology, Chinese Academy of Sciences, Shanghai 200050, China
University of Chinese Academy of Sciences, Beijing 100049, China  Siting Liu
School of Information Science and Technology, ShanghaiTech University,
Shanghai 201210, China  Shanghai Engineering Research Center of Energy
Efficient and Custom AI IC, ShanghaiTech University, Shanghai 201210, China
Vassilios Kovanis  Bradley Department of Electrical and Computer Engineering,
Virginia Tech Research Center, Arlington, Virginia 22203, USA  Cheng Wang
School of Information Science and Technology, ShanghaiTech University,
Shanghai 201210, China  Shanghai Engineering Research Center of Energy
Efficient and Custom AI IC, ShanghaiTech University, Shanghai 201210, China

###### Abstract

It has been found that noise-induced excitability in quantum well and quantum
dot semiconductor laser systems usually produce spike patterns of non-uniform
amplitude. In this letter, we experimentally record that an inter-subband
quantum cascade laser injected with a monochromatic laser exhibits a series of
highly-uniform spike trains in the time domain. Theoretical analysis
demonstrates that such high uniformity has its origin in the ultrashort
carrier lifetime of the quantum cascade laser gain medium that is typically
close to one picosecond.

††preprint: APS/123-QED

Excitability is a common phenomenon in various nonlinear dynamical systems,
including biological neurons, chemical reactions, and optical systems [1, 2].
When an excitable system is stimulated by a superthreshold perturbation, one
or several spikes are created [3]. This process is followed by a refractory
period, before which the system can not be excited again [4]. In the phase
space, the excitable system undergoes a large excursion from the stable
equilibrium [5]. Once the spike is fired, the system settles back to the
equilibrium state. On the other hand, a subthreshold perturbation can only
trigger a small response without producing any giant spikes. Semiconductor
lasers have been found to show excitability, when subjected to an external
control including saturable absorption [6], optical feedback [7], and the most
used optical injection [8, 9]. Optical injection unidirectionally injects the
monochromatic beam of light of a master laser into a slave laser with a
similar lasing frequency [10, 11, 12]. The slave diode laser is synchronized
with the master laser within the stable regime, which is bounded by the Hopf
bifurcation and the saddle-node bifurcation [13]. Spiking pulses are usually
observed in the vicinity of the saddle-node bifurcation, due to the existence
of homoclinic bifurcation [3]. The stimulus of noise randomly perturbs the
excitable laser systems and usually triggers a train of spikes [14, 15]. The
amplitudes of these noise-induced spikes are usually irregular and vary a lot
from one to another. On the other hand, finely controlled excitation of spikes
is valuable for developing spike-based neuromorphic computing systems, which
are attracting a lot of interest in recent years [16, 17, 18, 19].

Investigations of excitability are mostly based on near-infrared interband
semiconductor lasers, where the carrier lifetime is on the sub-nanosecond
scale. In contrast, the laser emission of mid-infrared and terahertz quantum
cascade lasers (QCLs) relies on the inter-subband transition. The carrier
lifetime of QCLs is around one picosecond, which is two to three orders of
magnitude smaller than common interband lasers [20]. Due to the ultrashort
carrier lifetime, QCLs have shown high stability with normal optical feedback,
while complex nonlinear dynamics became rare events [21, 22]. Nevertheless,
our previous work found that optical feedback with a tilted angle could
destabilize QCLs and trigger periodic oscillations, aperiodic oscillations,
and low-frequency oscillations [23, 24]. In addition, current modulation of
QCLs with optical feedback can stimulate low-frequency fluctuations and
extreme pulses [25, 26]. When subject to optical injection, previous reports
have theoretically identified both the Hopf bifurcation and the saddle-node
bifurcation of QCLs [27, 28]. Outside the stable locking regime, our recent
work demonstrated that QCLs mostly produced periodic oscillations [29]. In
this work, we show that a QCL is excitable in the vicinity of the saddle-node
bifurcation. It is strikingly found that the noise-induced spikes in the QCL
exhibit highly-uniform amplitude, which is in contrast to those in interband
semiconductor lasers. Theoretical modeling proves that the high uniformity of
the spike amplitudes inherently originates from the ultrashort carrier
lifetime of QCLs.

Figure 1: Experimental tabletop setup of the QCL injected optically with
monochromatic radiation. BS: beam splitter; FTIR: Fourier transform infrared
spectrometer, that has a resolution of 0.08 /cm; PD: photodetector; OSC: ultra
oscilloscope.

Figure 1 depicts the experimental setup for the QCL subject to the
monochromatic optical injection. The slave laser is a commercial distributed
feedback QCL (Thorlabs), that is driven by a continuous-wave electronic
current source. The operation temperature is maintained at 20 °C by using a
thermoelectric cooler. The master laser is a tunable external-cavity QCL
(Daylight solutions). The unidirectional injection is achieved through a
polarization-dependent isolator. The injection strength is controlled by a
half-wave plate and is monitored by a power meter. The optical spectrum is
measured by a Fourier transform infrared spectrometer (FTIR, Bruker) with a
resolution of 0.08 /cm. The optical signal is converted to the electrical one
through a HgCdTe photodetector (PD, Vigo) with a detection bandwidth of 560
MHz. The temporal waveform is recorded on a digital oscilloscope (OSC, 59 GHz
bandwidth), and the sampling rate is fixed at 5.0 GSample/s. The free-running
slave laser exhibits a lasing threshold of $I_{th}$=385 mA with an emission
wavenumber around 2182 /cm. The optical frequency detuning between the master
laser and the slave laser is finely tuned by adjusting the pump current of the
slave laser, instead of the master laser. The frequency tunability of the
slave laser is measured to be about -770.5 MHz/mA.

Figure 2: (a) Spike trains of the QCL recorded for several injection ratios.
(b) Mean amplitude and relative standard deviation of the spikes as a function
of the injection ratio.

When the slave QCL is pumped at 445 mA, the output power is 24.1 mW. The
optical detuning frequency is -2.45 GHz, which is located in the vicinity of a
saddle-node bifurcation. For injection ratios larger than -4 dB, the laser is
locked in the stable, phase-locked region and emits a continuous wave. When
the injection ratio is reduced to the range of -4 to -11 dB, the QCL produces
spikes as illustrated in Fig. 2(a). The spikes appear with random intervals,
due to the random nature of noise excitation. Surprisingly, the spike
amplitudes are highly uniform. This behavior is different from the noise-
induced spikes of common interband semiconductor lasers, where the amplitudes
vary substantially and distribute in a wide range [8]. Figure 2(a) also shows
that the occurrence rate of spikes depends strongly on the injection strength.
An injection ratio of -8 dB triggers the most number of spikes compared to
other injection powers. Figure 2(b) shows that the mean amplitude of the
spikes recorded in a duration of 1 ms rises nonlinearly with increasing
injection strength. Meanwhile, the relative standard deviation declines from
0.017 at the injection ratio of $R$=-11 dB down to 0.013 at $R$=-5 dB.
Therefore, strong optical injection improves the uniformity of the spikes.

Figure 3: Distribution of spikes for (a) intervals $>$ 0.1 $\mu\mathrm{s}$ and
(b) intervals $<$ 0.1 $\mu\mathrm{s}$. (c) Waveforms counted in (a), and (d)
waveforms counted in (b).

Figure 3 shows the distribution histogram of the spike intervals of the QCL,
where the injection ratio is -6 dB and the analyzed time duration is 1 ms. For
intervals above 0.1 $\mu\mathrm{s}$, Fig. 3(a) shows that the distribution
roughly follows the trend of exponential decay, which is a typical signature
of noise-induced excitability [14, 30]. The characteristic time of the decay
is about 1.3 $\mu\mathrm{s}$, which is governed by the noise strength. This
characteristic time is more than one order of magnitude larger than that of
common interband semiconductor lasers [14]. Therefore, kicking QCLs out of the
stable potential minimum is more difficult than the interband lasers. This
behavior happens because QCLs are intrinsically more stable than the latter,
owing to the ultrashort carrier lifetime and the small linewidth broadening
factor. Indeed, the high stability of QCLs has been widely demonstrated via
external perturbations from optical feedback [21, 22]. Temporal waveform
analysis in Fig. 3(c) shows that these intervals include three categories:
spike-to-spike interval (0.34 $\mu\mathrm{s}$), spike-to-burst interval (2.11
$\mu\mathrm{s}$ and 0.14 $\mu\mathrm{s}$), and burst-to-burst interval (1.24
$\mu\mathrm{s}$), which are illustrated in Fig. 3(c). For intervals below 0.1
$\mu\mathrm{s}$, Fig. 3(b) proves that the distribution roughly follows a
Gaussian function, with a mean value of 56 ns and a standard deviation of 4.7
ns. Temporal waveform analysis in Fig. 3(d) shows that this distribution comes
from the intra-burst pulses, which are evenly spaced. Each burst consists of a
different number of pulses, ranging from two up to more than ten.

In order to dissect the physical mechanism for producing spikes of uniform
amplitudes, we model the optically injected QCLs using a set of single-mode
rate equations, and the effects of quantum noise on the photon and carrier
dynamics are included as well [31]:

$\displaystyle\frac{dN_{3}}{dt}\\!=\\!\eta\frac{I}{q}\\!-\\!\frac{N_{3}}{\tau_{32}}\\!-\\!\frac{N_{3}}{\tau_{31}}\\!-\\!G_{0}S\Delta
N\\!+\\!F_{3}(t)$ (1)
$\displaystyle\frac{dN_{2}}{dt}\\!=\\!\frac{N_{3}}{\tau_{32}}\\!-\\!\frac{N_{2}}{\tau_{21}}\\!+\\!G_{0}S\Delta
N\\!+\\!F_{2}(t)$ (2)
$\displaystyle\frac{dN_{1}}{dt}\\!=\\!\frac{N_{3}}{\tau_{31}}\\!+\\!\frac{N_{2}}{\tau_{21}}\\!-\\!\frac{N_{1}}{\tau_{\text{out
}}}\\!+\\!F_{1}(t)$ (3)
$\displaystyle\frac{dS}{dt}\\!=\\!\left(\\!mG_{0}\Delta
N\\!-\\!\frac{1}{\tau_{p}}\\!\right)\\!S\\!+\\!m\beta\frac{N_{3}}{\tau_{sp}}\\!+\\!2k_{c}\sqrt{S_{inj}S}\cos\phi+F_{s}(t)$
(4)
$\displaystyle\frac{d\phi}{dt}\\!=\\!\frac{\alpha_{H}}{2}\\!\left(\\!mG_{0}\Delta
N\\!-\\!\frac{1}{\tau_{p}}\\!\right)\\!-\\!2\pi\Delta
F_{inj}\\!-\\!k_{c}\sqrt{\frac{S_{inj}}{S}}\sin\phi\\!+\\!F_{\phi}(t)$ (5)

where $N_{3,2,1}$ are carrier numbers in the upper lasing level, the lower
lasing level, and the bottom level, respectively. $S$ is the photon number in
the injection-locked laser cavity, and $\phi$ is the phase difference between
the slave laser and the master laser. $I$ is the pump current, $\eta$ is the
current injection efficiency, $G_{0}$ is the gain coefficient
($G_{0}\\!=\\!5.3\\!\times\\!10^{4}$ /s), $m$ is the number of gain stage
($m$=30), $\beta$ is the spontaneous emission factor ($\beta\\!=\\!10^{-6}$),
and $\Delta N$ is the population inversion expressed as $\Delta
N\\!=\\!N_{3}\\!-\\!N_{2}$. $S_{inj}$ is the photon number of the master
laser, and the injection ratio is defined as $R_{inj}\\!=\\!S_{inj}/S_{0}$,
with $S_{0}$ being the photon number of the free-running slave laser. $\Delta
F_{inj}$ is the optical detuning frequency between the master laser and the
slave laser, and $k_{c}$ is the coupling coefficient between the two lasers.
The carrier lifetimes for each level are $\tau_{32}$=2.0 ps, $\tau_{31}$=2.4
ps, $\tau_{21}$=0.5 ps, $\tau_{out}$=0.5 ps. The photon lifetime is
$\tau_{p}$=3.7 ps, and the spontaneous emission lifetime is 7.0 ns. The
carrier noise $F_{3,2,1}$ and the spontaneous emission noise, quantum noise
$F_{S,\phi}$ are modeled by the Langevin noise sources, respectively [32].
Detailed correlation relations between the noise sources refer to [31]. The
linewidth broadening factor is fixed at $\alpha_{H}$=-0.5. We note that
measured linewidth broadening factors of QCLs in the literature range from
-1.0 up to more than 2.0 [21, 33]. However, the optical injection can reduce
its value due to the reduction of carrier population [34].

Figure 4: Simulated spike trains of (a) the common QCL with ultrashort carrier
lifetimes and (b) the artificial QCL with 100-time longer carrier lifetimes.
The detuning frequency is -2.1 GHz in (a) and -1.2 GHz in (b), respectively.
The injection ratio is -5 dB for both cases.

In the vicinity of the saddle-node bifurcation, the simulation in Fig. 4(a)
successfully reproduces the spike trains with highly uniform amplitudes, which
agrees well with the experimental results in Fig. 2(a). The relative standard
deviation of the simulated spikes is 0.007. This deviation is smaller than the
experimental one in Fig. 2(b). It is likely due to the fact that the rate
equation model does not take into account technical noise sources, including
the current source noise, the operation temperature fluctuation, and the
mechanical vibration. Once we remove the noise sources in the rate equation
model, the spiking dynamics do not occur. Therefore, it is verified that the
generation of spikes in QCLs arises from the noise-induced excitability. In
order to explore the mechanism governing the uniformity of amplitudes, we
artificially enlarge all the carrier lifetimes by a factor of 100, while other
parameters remain unchanged. That is, the carrier lifetime of the QCL becomes
about 100 ps. Figure 4(b) shows the spike trains of this artificial QCL with
optical injection. Interestingly, the spike amplitudes are not uniform anymore
but vary significantly from spike to spike. The corresponding relative
standard deviation reaches as high as 0.07, which is one order of magnitude
larger than the one with ultrashort carrier lifetimes in Fig. 4(a). These non-
uniform spikes are similar to those commonly observed in interband
semiconductor lasers with a lifetime of sub-nanosecond [8]. Therefore, we can
conclude that the high uniformity of the spikes in QCLs is owing to the
ultrashort carrier lifetime. Indeed, the fast gain recovery process in the QCL
gain medium is helpful to store strong enough energy for each firing spike,
whenever the noise stimulus perturbs the laser system [35, 36].

We remark that the uniform spikes occur in a broad regime (bounded by the
injection ratio and the detuning frequency) in the vicinity of the saddle-node
bifurcation. In addition, we have observed uniform spikes for the pump
currents both close to and well above the lasing threshold of the QCL.
Therefore, the excitation of uniform spikes in QCLs with optical injection is
robust and reliable. These uniform spikes are highly valuable for developing
neuromorphic computing systems or spiking neural networks. These networks
consist of many spiking neurons and the neuron dynamics are usually described
by the popular leaky-integrate-and-fire (LIF) model [1]. In this model, the
membrane potential of the neuron continuously experiences an exponential decay
toward the resting value. Once the neuron receives a spike, the membrane
potential increases by an amount proportional to the amplitude of the spike.
When the potential reaches a certain threshold, the neuron fires a spike and
resets its membrane potential. With this repeated process, the neuron
generates a series of spikes or a spike train. To perform neuromorphic
computing using these neurons, the spike train is usually encoded by the
spiking rate (rate coding) or the spiking time (temporal coding) [37]. To
successfully code the message, it is crucial to ensure the uniformity of the
spikes, i.e., all the input spikes and output spikes must have identical
amplitude. Otherwise, a large input spike may incorrectly lead to the fast
rate or early appearance of the output spike, and vice versa. Consequently, an
error occurs due to the change of the message encoded on the output spike
train. Indeed, CMOS-based electronic neuromorphic computing systems usually
own uniform spike trains [38, 39]. However, the spikes produced from near-
infrared lasers like vertical-cavity surface-emitting lasers are usually not
uniform, which prevents high-quality implementation of the LIF model based on
laser diodes [19, 40].

In summary, we have experimentally recorded the observation of highly uniform
spike trains in a tabletop optically injected quantum cascade system. The
spike trains are excited via quantum noise, in the vicinity of a saddle-node
bifurcation, where a collision and disappearance of two equilibria points are
happening. Our theoretical analysis demonstrates that such high uniformity is
born from the ultra-short carrier lifetime of the optical gain medium, which
is more than two orders of magnitude smaller than that of common quantum well
laser diodes. We assert that uniform spikes are of prime importance for
developing next-generation photonic neuromorphic computing systems. Future
work will investigate the detailed characteristics of spikes excited by
properly controlled input signals instead of just quantum noise, such as the
excitation threshold, the spike latency, and the refractory period.

This work was financially supported by Shanghai Natural Science Foundation
(20ZR1436500) and the work of VK was supported via generous gifts to the
Virginia Tech Foundation.

## References

  * Izhikevich [2007] E. M. Izhikevich,  _Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting_ (MIT press, Cambridge, MA, 2007).
  * Izhikevich [2004] E. M. Izhikevich, Which model to use for cortical spiking neurons?, IEEE Trans. Neural Netw. 15, 1063 (2004).
  * Wieczorek  _et al._ [2002] S. Wieczorek, B. Krauskopf, and D. Lenstra, Multipulse excitability in a semiconductor laser with optical injection, Phys. Rev. Lett. 88, 063901 (2002).
  * Garbin  _et al._ [2017] B. Garbin, A. Dolcemascolo, F. Prati, J. Javaloyes, G. Tissoni, and S. Barland, Refractory period of an excitable semiconductor laser with optical injection, Phys. Rev. E 95, 012214 (2017).
  * Krauskopf  _et al._ [2003] B. Krauskopf, K. Schneider, J. Sieber, S. Wieczorek, and M. Wolfrum, Excitability and self-pulsations near homoclinic bifurcations in semiconductor laser systems, Opt. Commun. 215, 367 (2003).
  * Barbay  _et al._ [2011] S. Barbay, R. Kuszelewicz, and A. M. Yacomotti, Excitability in a semiconductor laser with saturable absorber, Opt. Lett. 36, 4476 (2011).
  * Giudici  _et al._ [1997] M. Giudici, C. Green, G. Giacomelli, U. Nespolo, and J. R. Tredicce, Andronov bifurcation and excitability in semiconductor lasers with optical feedback, Phys. Rev. E 55, 6414 (1997).
  * Goulding  _et al._ [2007] D. Goulding, S. P. Hegarty, O. Rasskazov, S. Melnik, M. Hartnett, G. Greene, J. G. McInerney, D. Rachinskii, and G. Huyet, Excitability in a quantum dot semiconductor laser with optical injection, Phys. Rev. Lett. 98, 153903 (2007).
  * Dillane  _et al._ [2019] M. Dillane, J. Robertson, M. Peters, A. Hurtado, and B. Kelleher, Neuromorphic dynamics with optically injected quantum dot lasers, Eur. Phys. J. B 92, 197 (2019).
  * Valagiannopoulos and Kovanis [2021] C. Valagiannopoulos and V. Kovanis, Injection-locked photonic oscillators: legacy results and future applications, IEEE Antennas Propag. Mag. 63, 51 (2021).
  * Himona  _et al._ [2022] G. Himona, V. Kovanis, and Y. Kominis, Isochrons, phase response and synchronization dynamics of tunable photonic oscillators, Phys. Rev. Res. 4, L012039 (2022).
  * Herrera  _et al._ [2022] D. J. Herrera, V. Kovanis, and L. F. Lester, Using transitional points in the optical injection locking behavior of a semiconductor laser to extract its dimensionless operating parameters, IEEE J. Sel. Top. Quantum Electron. 28, 1800109 (2022).
  * Ohtsubo [2017] J. Ohtsubo,  _Semiconductor Lasers: Stability, Instability and Chaos_ (Springer, Cham, Switzerland, 2017).
  * Kelleher  _et al._ [2011a] B. Kelleher, D. Goulding, G. Huyet, E. A. Viktorov, T. Erneux, and S. P. Hegarty, Dimensional signature on noise-induced excitable statistics in an optically injected semiconductor laser, Phys. Rev. E 84, 026208 (2011a).
  * Kelleher  _et al._ [2011b] B. Kelleher, C. Bonatto, G. Huyet, and S. P. Hegarty, Excitability in optically injected semiconductor lasers: contrasting quantum-well-and quantum-dot-based devices, Phys. Rev. E 83, 026207 (2011b).
  * Turconi  _et al._ [2013] M. Turconi, B. Garbin, M. Feyereisen, M. Giudici, and S. Barland, Control of excitable pulses in an injection-locked semiconductor laser, Phys. Rev. E 88, 022923 (2013).
  * Prucnal  _et al._ [2016] P. R. Prucnal, B. J. Shastri, T. F. de Lima, M. A. Nahmias, and A. N. Tait, Recent progress in semiconductor excitable lasers for photonic spike processing, Adv. Opt. Photonics 8, 228 (2016).
  * Brunner  _et al._ [2018] D. Brunner, B. Penkovsky, B. A. Marquez, M. Jacquot, I. Fischer, and L. Larger, Tutorial: photonic neural networks in delay systems, J. Appl. Phys. 124, 152004 (2018).
  * Robertson  _et al._ [2019] J. Robertson, E. Wade, Y. Kopp, J. Bueno, and A. Hurtado, Toward neuromorphic photonic networks of ultrafast spiking laser neurons, IEEE J. Sel. Top. Quantum Electron. 26, 7700715 (2019).
  * Faist  _et al._ [1994] J. Faist, F. Capasso, D. L. Sivco, C. Sirtori, A. L. Hutchinson, and A. Y. Cho, Quantum cascade laser, Science 264, 553 (1994).
  * Zhao  _et al._ [2020] B. B. Zhao, X. G. Wang, and C. Wang, Strong optical feedback stabilized quantum cascade laser, ACS Photonics 7, 1255 (2020).
  * Wang  _et al._ [2020] X. G. Wang, B. B. Zhao, F. Grillot, and C. Wang, Spectral linewidth reduction of quantum cascade lasers by strong optical feedback, J. Appl. Phys. 127, 073104 (2020).
  * Wang  _et al._ [2021] X. G. Wang, B. B. Zhao, Y. Deng, V. Kovanis, and C. Wang, Nonlinear dynamics of a quantum cascade laser with tilted optical feedback, Phys. Rev. A 103, 023528 (2021).
  * Zhao  _et al._ [2021] B. B. Zhao, X. G. Wang, and C. Wang, Low-frequency oscillations in quantum cascade lasers with tilted optical feedback, IEEE J. Sel. Top. Quantum Electron. 28, 1800207 (2021).
  * Jumpertz  _et al._ [2016] L. Jumpertz, K. Schires, M. Carras, M. Sciamanna, and F. Grillot, Chaotic light at mid-infrared wavelength, Light: Sci. Appl. 5, e16088 (2016).
  * Spitz  _et al._ [2020] O. Spitz, J. G. Wu, A. Herdt, G. Maisons, M. Carras, W. E. Elsäßer, C. W. Wong, and F. Grillot, Extreme events in quantum cascade lasers, Adv. Photonics 2, 066001 (2020).
  * Wang  _et al._ [2013] C. Wang, F. Grillot, V. Kovanis, J. D. Bodyfelt, and J. Even, Modulation properties of optically injection-locked quantum cascade lasers, Opt. Lett. 38, 1975 (2013).
  * Erneux  _et al._ [2013] T. Erneux, V. Kovanis, and A. Gavrielides, Nonlinear dynamics of an injected quantum cascade laser, Phys. Rev. E 88, 032907 (2013).
  * Peng  _et al._ [2022] Y. B. Peng, B. B. Zhao, and C. Wang, Nonlinear dynamics of a quantum cascade laser with optical injection, Opt. Express 30, 27593 (2022).
  * Kramers [1940] H. A. Kramers, Brownian motion in a field of force and the diffusion model of chemical reactions, Physica 7, 284 (1940).
  * Wang  _et al._ [2018] X. G. Wang, B. B. Zhao, F. Grillot, and C. Wang, Frequency noise suppression of optical injection-locked quantum cascade lasers, Opt. Express 26, 15167 (2018).
  * Coldren  _et al._ [2012] L. A. Coldren, S. W. Corzine, and M. L. Mashanovitch,  _Diode Lasers and Photonic Integrated Circuits_ (Wiley, Hoboken, NJ, 2012).
  * Opačak  _et al._ [2021] N. Opačak, F. Pilat, D. Kazakov, S. Dal Cin, G. Ramer, B. Lendl, F. Capasso, and B. Schwarz, Spectrally resolved linewidth enhancement factor of a semiconductor frequency comb, Optica 8, 1227 (2021).
  * Wang  _et al._ [2015] C. Wang, M. E. Chaibi, H. M. Huang, D. Erasme, P. Poole, J. Even, and F. Grillot, Frequency-dependent linewidth enhancement factor of optical injection-locked quantum dot/dash lasers, Opt. Express 23, 21761 (2015).
  * Choi  _et al._ [2008] H. Choi, L. Diehl, Z. K. Wu, M. Giovannini, J. Faist, F. Capasso, and T. B. Norris, Gain recovery dynamics and photon-driven transport in quantum cascade lasers, Phys. Rev. Lett. 100, 167401 (2008).
  * Green  _et al._ [2009] R. P. Green, A. Tredicucci, N. Q. Vinh, B. Murdin, C. Pidgeon, H. E. Beere, and D. A. Ritchie, Gain recovery dynamics of a terahertz quantum cascade laser, Phys. Rev. B 80, 075303 (2009).
  * Gerstner and Kistler [2002] W. Gerstner and W. M. Kistler,  _Spiking Neuron Models: Single Neurons, Populations, Plasticity_ (Cambridge University Press, Cambridge, UK, 2002).
  * El Srouji  _et al._ [2022] L. El Srouji, A. Krishnan, R. Ravichandran, Y. Lee, M. On, X. Xiao, and S. J. Ben Yoo, Photonic and optoelectronic neuromorphic computing, APL Photonics 7, 051101 (2022).
  * Chalkiadakis and Hizanidis [2022] D. Chalkiadakis and J. Hizanidis, Dynamical properties of neuromorphic Josephson junctions, Phys. Rev. E 106, 044206 (2022).
  * Owen-Newns  _et al._ [2023] D. Owen-Newns, J. Robertson, M. Hejda, and A. Hurtado, GHz rate neuromorphic photonic spiking neural network with a single vertical-cavity surface-emitting laser (VCSEL), IEEE J. Sel. Top. Quantum Electron. 29, 1500110 (2023).