IISER]Department of Physics, Indian Institute of Science Education and
Research (IISER), Pune 411008, India IISER]Department of Physics, Indian
Institute of Science Education and Research (IISER), Pune 411008, India

# Simultaneous detection of spin and orbital angular momentum of light through
scattering from a single silver nanowire

Diptabrata Paul<EMAIL_ADDRESS>[  Deepak K. Sharma
Laboratoire Interdisciplinaire Carnot de Bourgogne, UMR 6303 CNRS, Université
de Bourgogne Franche-Comté, 9 Avenue Alain Savary, 21000 Dijon, France  G V
Pavan Kumar<EMAIL_ADDRESS>[

###### Abstract

In recent times the spin angular momentum (SAM) and orbital angular momentum
(OAM) of light have gained prominence because of their significance in optical
communication systems, micromanipulation, sub-wavelength position sensing. To
this end, simultaneous detection of SAM and OAM of light beam is one of the
important topics of research from both application and fundamental spin-orbit
interaction (SOI) point of view. While interferometry and metasurface based
approaches have been able to detect the states, our approach involves elastic
scattering from a monocrystalline silver nanowire for the simultaneous
detection of SAM and OAM state of a circularly polarized Laguerre-Gaussian
(LG) beam. By employing Fourier plane (FP) microscopy, the transmitted
scattered light intensity distribution in the FP is analyzed to reconstruct
the SAM and OAM state unambiguously. The SAM and OAM induced transverse energy
flow as well as the polarization dependent scattering characteristics of the
nanowire is investigated to understand the detection mechanism. Our method is
devoid of complex nanofabrication techniques required for metasurface based
approaches and to our knowledge, is a first example of single nano-object
based simultaneous SAM and OAM detection. The study will further the
understanding of SOI effects and can be useful for on-chip optical detection
and manipulation.

##  1 Introduction

A light beam can possess spin angular momentum (SAM) manifested through its
circular polarization state, as well as orbital angular momentum (OAM),
manifested through its helical phase front 1, 2, 3. SAM can have two
orthogonal states corresponding to left- and right-circularly polarized light,
whereas OAM is unbounded with its infinite orthogonal states, corresponding to
the topological charge of vortex beams such as Laguerre-Gaussian beams
($\textrm{LG}^{n}_{m}$, $n$ topological charge, $m$ radial index) 4, 5. This
angular momentum of light corresponds to an energy flow in the beam in the
transverse plane with respect to the propagation direction, termed as
transverse energy flow 6. Interaction of the transverse energy flow of the
beam with a medium may lead to spin-orbit interaction (SOI) induced effects 7,
8, 9, 10 such as directional scattering 11, interconversion of angular
momentum 12, 13 as well as micro manipulation 14 and position sensing 15.
Additionally, beams possessing SAM and OAM have been envisaged for optical
communication 16, 17, 18, 19, quantum information 20 applications due to their
orthogonal states.

To this end, deterministic sensing of SAM and OAM has become one of the
important topics of research in recent years. Detection of OAM in macroscopic
scale has been achieved through conventional interferometry 21, 22, 23 as well
as diffraction 24, 25, 26, 27 and projection through a reciprocal forked
diffraction grating 20. These methods are constrained by their complicated
optical setups. With the advancement of nanofabrication, generation as well as
detection of OAM and SAM has been achieved at nanoscale 28, 29, 30, 31, 32. At
these length scales, the on chip plasmonic devices offer great advancement
from the interferometric based techniques with their lower footprint. But, in
most cases these devices involve complex fabrication methods.

Figure 1: Schematic of the experimental configuration. An AgNW placed on a
glass substrate scatters the incident LG beam to determine its total angular
momentum state.

Thus, there is a requirement of methods which are devoid of these
complications and possess smaller footprint for chip-scale applications.
Additionally, differential sensing of spin and orbital angular momentum
component of a light beam using a single nano-object is yet to be achieved.
With this motivation, herein we propose a chemically synthesized quasi-one-
dimensional monocrystalline silver nanowire (AgNW) as a sensor for the
simultaneous detection of SAM and OAM of light. AgNWs have gained prominence
because of their simple preparation method and monocrystalline
characteristics. The solution phase synthesis allows them to be self-assembled
on any desired substrate and having necessary tunability 33, 34. The uniform
morphology renders them useful for low loss surface plasmon waveguiding
property 35 and makes them ideal for fundamental photonic circuit applications
such as plasmonic waveguide 35, 36, plasmonic antenna 37, 38, logic gates 39
as well as for opto-electrical applications 40. Additionally, in our previous
reports we have shown how an AgNW can act as a plasmonic scatterer in sensing
total energy flow of an incident beam 41, 42 and as well as probing the
inherent spin-Hall effect 43.

In this work, a single AgNW dropcasted on a glass substrate has been used as a
nanoscopic strip diffractor and the scattered pattern has been extensively
analyzed to ascertain the spin and orbital angular momentum state of the
incident light, as shown in Figure 1. By considering the total transverse
energy flow, as well as the individual component due to spin flow and orbital
flow at the focus, we are able to distinguish their individual effect on the
scattering pattern, allowing us to unambiguously determine their value.

##  2 Results and Discussion

###  2.1 Theoretical Calculation

In order to distinguish spin and orbital angular momentum, we need to
characterize their contribution to the energy flow of the beam 6. The energy
flow of an electromagnetic field, defined by its Poynting vector, can be
decomposed into its longitudinal and transverse part. For example, Laguerre-
Gaussian (LG) beams, commonly known as vortex beams, possess transverse energy
flow corresponding to its orbital flow, owing to its helical wave front. The
magnitude and the sense of flow gets dictated by the magnitude and sign of the
topological charge of the beam. On the other hand, circularly polarized
Gaussian beams exhibit transverse energy flow corresponding to its spin flow,
the sense of which gets dictated by the handedness of the circular
polarization.

Figure 2: Theoretically calculated focal orbital flow density (OFD), spin flow
density (SFD) and total transverse energy flow distribution for (a) left-
circularly polarized ($\sigma_{+}$) and, (b) right-circularly polarized
($\sigma_{-}$) $\textrm{LG}^{1}_{0}$ beam. Inset shows the corresponding total
field intensity. (c) and (d) represents the OFD, SFD and total transverse
energy flow distribution line profile along $x$ axis from $0$ to $2\lambda$.

To investigate the interaction due to transverse energy flow at nanoscale, we
employ the Debye-Wolf integral formulation 44, 45 to find out the optical
electric field ($\mathbf{E}=(E_{x},E_{y},E_{z})$) and magnetic field
($\mathbf{H}=(H_{x},H_{y},H_{z})$) components at the focal plane of a $0.5$ NA
lens 46. From these, the corresponding spin flow density (SFD,
$\mathbf{{p}^{\perp}_{s}}$) and orbital flow density (OFD,
$\mathbf{{p}^{\perp}_{o}}$) about $z$ axis (propagation direction) can be
calculated by considering the transverse component (denoted by superscript
$\perp$) of the following equations 6:

$\displaystyle\mathbf{p_{s}}=\frac{g}{4\omega}\textrm{Im}[\nabla\times(\varepsilon_{0}(\mathbf{E^{*}}\times\mathbf{E}))+\mu_{0}(\mathbf{H^{*}}\times\mathbf{H}))],$
(1)

$\displaystyle\mathbf{p_{o}}=\frac{g}{2\omega}\textrm{Im}[\varepsilon_{0}((\mathbf{E^{*}}\cdot(\nabla)\mathbf{E})+\mu_{0}(\mathbf{H^{*}}\cdot(\nabla)\mathbf{H})].$
(2)

Here $g=(8\pi)^{-1}$ and $\omega$ is the oscillating frequency,
$\varepsilon_{0}$ and $\mu_{0}$ are the electric permittivity and magnetic
permeability of vaccum respectively and
$\mathrm{{p}^{\perp}_{o,i}}\propto\textrm{Im}[\Sigma_{j}{E}^{*}_{j}(\nabla_{i}){E}_{j}]$
($i,j=(x,y,z)$). Figure 2(a) and (b) show the calculated OFD, SFD and total
transverse energy flow density
($\mathbf{p}^{\perp}=\mathbf{p_{s}}^{\perp}+\mathbf{p_{o}}^{\perp}$)
corresponding to left-circularly polarized ($\sigma_{+}$) and right-circularly
polarized ($\sigma_{-}$) $\textrm{LG}^{1}_{0}$ beams at wavelength
$\lambda=633$ nm. The total electric field intensity at the focal plane is
shown in the inset. As expected, for both $\sigma_{+}$ and $\sigma_{-}$
polarized $\textrm{LG}^{1}_{0}$ beam, the OFD has similar magnitude and
similar counter-clock orbital flow sense. But, for $\sigma_{+}$
$\textrm{LG}^{1}_{0}$ beam, the SFD and the corresponding spin flow have both
clockwise and counter-clock energy flow in its inner and outer part
respectively. Additionally, the magnitude of SFD is higher in its inner radius
($|\mathbf{p_{s}}^{\perp}|=0.1085$ at $R=0.322\lambda$) and slightly lower in
its outer radius ($|\mathbf{p_{s}}^{\perp}|=0.0731$ at $R=1.068\lambda$). It
is exactly opposite for the $\sigma_{-}$ $\textrm{LG}^{1}_{0}$ beam. Hence,
for $\sigma_{+}$ $\textrm{LG}^{1}_{0}$ beam, counter-clock orbital flow
towards the center cancels the clockwise spin flow, whereas, in the outer part
the counter-clock spin flow gets added up with orbital flow, resulting an
overall counter clock total transverse energy flow having a maximum magnitude
of $|\mathbf{p}^{\perp}|=0.1515$ at radius $R=0.868\lambda$. For $\sigma_{-}$
$\textrm{LG}^{1}_{0}$ the counter-clock inner spin flow gets added up with the
counter-clock orbital flow resulting in higher counter-clock transverse energy
flow, $|\mathbf{p}^{\perp}|=0.2551$ at $R=0.4\lambda$. The OFD, SFD and the
total transverse energy flow density radial line profile at the focal plane
from $0$ to $2\lambda$ corresponding to $\sigma_{+}$ and $\sigma_{-}$
polarized $\textrm{LG}^{1}_{0}$ beams are shown in figure 2(c) and (d)
respectively. In this context, it is important to point out the subtle
differences of the intensity profile for $\sigma_{+}$ and $\sigma_{-}$
polarized $\textrm{LG}^{1}_{0}$ beams as focusing of circularly polarized LG
beam aids spin-to-orbital angular momentum transformation 12. The field
calculations and corresponding phase profiles are shown in the supplementary
information S1.

###  2.2 Experimental Implementation

The focal optical field and its orbital and spin flow constituents, upon
interaction with a nano structure will lead to flow dependant scattering
allowing us to determine its state. Although different nano-structures have
been shown to have exhibited transverse energy flow dependant scattering 30,
47, 31, 41, 42, differential sensing of spin and orbital flow is yet to be
achieved by a single nano-object. To that end, we use a single crystalline
silver nanowire (AgNW) to scatter the incident optical field and sense its
properties. The implementation of forward scattering configuration is achieved
by dropcasting AgNWs of width $\approx 350$ nm onto a glass cover-slip, which
acts as a nanoscopic strip scatterer 24, 27 for the incident optical fields. A
$50\times 0.5$ NA objective lens has been used to focus the incident beam and
the scattered signal was collected using a $100\times 1.49$ NA lens, in the
transmission configuration (Figure 3(a)). A low power of $0.13\mathrm{\mu W}$
has been used for the scattering of $\textrm{LG}^{\pm 1}_{0}$ beam. Fourier
plane (FP) intensity distribution of the forward scattered light was collected
by projecting the back focal plane of the collection objective lens onto a CCD
using relay optics 48, 49. The collected FP intensity distribution can be
divided into two regions, sub-critical region ($\textrm{NA}<1$) and super-
critical region ($\textrm{NA}>1$). The incident NA part (NA $0-0.5$) in the
sub-critical region is dominated by the un-scattered incident light, and is
omitted from our analysis, appears as a black disk at the center of the
measured FP. Details of the experimental configuration is elaborated in the
experimental section and supplementary information S2.

Figure 3: (a) Experimental configuration of measurement scheme for the Fourier
plane (FP). Inset shows the bright field optical micrograph of AgNW.
Experimentally measured Fourier plane intensity distribution for scattering of
(b) $x$ polarized and, (c) $y$ polarized $\textrm{LG}^{1}_{0}$ beam. Inner and
outer white circles represent the critical angle at the air-glass interface
and the collection limit of the objective lens NA respectively. Incident NA
region is blocked and appears as a black disk at the center. (d) Simulated
near-field electric field intensity distribution at the air-glass interface
for $x$ and $y$ polarized incident beam. (e) Corresponding scattered total
transverse energy flow density ($\textrm{p}^{\perp}_{scat}$) distribution for
$x$ and $y$ polarized $\textrm{LG}^{1}_{0}$ beam. The black arrows indicate
the energy flow direction.

###  2.3 Polarization dependence of scattering intensity

To determine the contribution due SAM and OAM in the scattered FP, it is
imperative to understand the scattered FP patterns of its constituent linear
poalrization components. The longitudinally polarized optical field ($E_{x}$)
with respect to the nanowire long axis ($x$ axis) gets scattered in the sub-
critical region, as shown in the experimentally measured data in figure 3(b)
for $x$ polarized $\textrm{LG}^{1}_{0}$ beam. On the other hand, the
transversely polarized fields ($E_{y}$) lead to near-field (NF) at the NW
edges, and consequent signal at the super-critical region in the FP as well as
sub-critical region, shown in figure 3(c) for $y$ polarized
$\textrm{LG}^{1}_{0}$ beam. Additionally, the polarization of the scattered
light for the $x$ polarized incident beam is predominantly $x$ polarized and
similarly, for $y$ polarized beam scattering, the light scattered in the sub-
and super-critical region is mostly $y$ polarized (supplementary information
S3). Thus, for scattering of circularly polarized $\textrm{LG}^{1}_{0}$ beam,
the state of circular polarization can be determined by analyzing the sub-
critical region intensity distribution as only in this regime both in-plane
$x$ and $y$ components of polarization scatter. Whereas, in super-critical
region, predominantly scattered light is from $y$ polarized field and will
posses information about OAM state of the incident beam. This has been further
confirmed by analyzing the FP due to $\sigma_{+}$ polarized
$\textrm{LG}^{1}_{0}$ beam (supplementary information S3), where upon
analyzing the FP intensity distribution for $x$ and $y$ polarized components,
we can recover the FP intensity distribution for $x$ and $y$ polarized
intensity distribution as shown in figure 3(b) and (c) respectively. It must
also be observed that both FP intensity distribution for $x$ and $y$ polarized
$\textrm{LG}^{1}_{0}$ beam exhibit two lobes as well as preferential
scattering pattern due to their orbital flow 42. But, the preference for $x$
polarized case is opposite to that of $y$ polarized case. This can be
attributed to the polarization dependent interaction of the NW with the
incident beam as shown through numerical simulation in figure 3(d) and (e).
The details of the finite element method numerical simulation and geometry is
given in supplementary information S4. For $x$ polarized $\textrm{LG}^{1}_{0}$
beam the field gets scattered from the NW and the coresponding total
transverse energy flow ($\textrm{p}^{\perp}_{scat}$) points away from the NW
and opposite to the incident beam $\textrm{p}^{\perp}$. On the other hand, $y$
polarized $\textrm{LG}^{1}_{0}$ beam aids the interaction of the transverse
energy flow with the nanowire and leads to generation of evanescent near-field
at nanowire edges, apparent with the presence of scattered light signal in the
super-critical region of FP in figure 3(c). Consequently, the scattered field
$\textrm{p}^{\perp}_{scat}$ flow is along the same direction as that of the
incident beam. Thus, a nanowire through its scattering can detect OAM modes as
well as its linear polarization states. The effect is qualitatively similar to
a nanoscopic strip scatterer 24, 27.

Figure 4: Experimentally measured scattered FP intensity distribution of (a)
$\sigma_{+}$ polarized and (b) $\sigma_{-}$ polarized $\textrm{LG}^{1}_{0}$
beam. The corresponding intensity distribution profile in the sub-critical
region of the $k_{y}>0$ part of the FP for $\sigma_{+}$ and $\sigma_{-}$
polarized $\textrm{LG}^{1}_{0}$ beam is shown in (c) and (d) respectively. The
same for super-critical region is shown in (e) and (f). The black arrows
indicate the effective directionality.

###  2.4 Nanowire as a sensor for OAM and SAM state

Next, we investigate the FP intensity distribution for $\sigma_{+}$ and
$\sigma_{-}$ polarized $\textrm{LG}^{\pm 1}_{0}$ beam to determine OAM and SAM
state of the beam. Figure 4(a) and (b) shows the FP intensity distribution for
$\sigma_{+}$ and $\sigma_{-}$ $\textrm{LG}^{1}_{0}$ beam respectively. We
analyze the pattern by dividing the FP into four parts and looking into the
intensity distribution of either upper ($k_{y}>0$ region) or lower half
($k_{y}<0$ region). Additionally, the SAM and OAM signature can be obtained by
analysing the signal in the sub-critical and super-critical region intensity
distribution respectively. Figure 4(c) and (d) shows the intensity
distribution profile within the sub-critical region of the $k_{y}>0$ part in
the FP, for $\sigma_{+}$ and $\sigma_{-}$ polarized $\textrm{LG}^{1}_{0}$ beam
respectively. The preference of the scattering is opposite for $\sigma_{+}$
polarized $\textrm{LG}^{1}_{0}$ beam scattering with respect to that of
$\sigma_{-}$ $\textrm{LG}^{1}_{0}$ beam data as indicated by the black arrows.
The preference is the result of spin flow. Since the dominant contribution tot
he spin flow occurs due to the inner circular part of the beam as shown in
figure 2, the preference follows the same handedness. The super-critical
region intensity distribution profile in the $k_{y}>0$ half for both
$\sigma_{+}$ and $\sigma_{-}$ polarized $\textrm{LG}^{1}_{0}$ beam exhibit
similar preference as shown in Figure 4 (e) and (f). This preference as well
as the number of lobes indicates the OAM state of the incident beam 42.

Figure 5: Experimentally measured scattered FP intensity distribution (a)
$\sigma_{+}$ and (b) $\sigma_{-}$ $\textrm{LG}^{-1}_{0}$ beam. The
corresponding intensity distribution profile in the sub-critical region of the
$k_{y}>0$ part of the FP for $\sigma_{+}$ and $\sigma_{-}$ polarized
$\textrm{LG}^{-1}_{0}$ beam is shown in (c) and (d) respectively. The same for
super-critical region is shown in (e) and (f). The black arrows indicate the
effective directionality.

The intensity distributions can be further quantified by calculating the
directionality, $D=(I_{2}-I_{1})/(I_{2}+I_{1})$, where $I_{i}$ represents the
integrated intensity of $i^{th}$ lobe for sub- or super-critical region. From
the experimentally measured data in figure 4 (c) and (d) we find,
$D_{sub}=0.094$ for $\sigma_{+}$ $\textrm{LG}^{1}_{0}$ beam and
$D_{sub}=-0.123$ for $\sigma_{-}$ $\textrm{LG}^{1}_{0}$ beam in the sub-
critical region. Thus, $D_{sub}>0$ for $\sigma_{+}$ polarization case and
$D_{sub}<0$ for $\sigma_{-}$ polarized beam, i.e., sign of the directionality
in the sub-critical region inverts depending upon the spin flow direction and
hence the SAM state. On the other hand, the preferential scattering in the
super-critical region of the FP indicates the sense of orbital flow and hence
the OAM state: $D_{super}=-0.284$ for $\sigma_{+}$ polarized and
$D_{super}=-0.563$ for $\sigma_{-}$ polarized $\textrm{LG}^{1}_{0}$ beam,
i.e., directionality sign ($D_{super}<0$) remains the same. Hence, $D_{sub}>0$
and $D_{super}<0$ indicates $\sigma_{+}$ $\textrm{LG}^{1}_{0}$ beam excitation
and $D_{sub}<0$ and $D_{super}<0$ indicates $\sigma_{-}$ $\textrm{LG}^{1}_{0}$
beam excitation. Although the differential sensing of the orbital and spin
flow is achieved, the total transverse energy flow influences the effective
directionality magnitude in the super-critical region to some extent. The
measured directionality for the super-critical region is significantly lower
for $\sigma_{+}$ $\textrm{LG}^{1}_{0}$ beam than for $\sigma_{-}$
$\textrm{LG}^{1}_{0}$ beam, which is due to the higher total transverse energy
flow density for $\sigma_{-}$ polarized beam, as shown in figure 2. For sub-
critical region the change is much smaller due lower directionality value as
well as positional sensitivity. The directionality values will be exactly
opposite in sign when we consider the lower half ($k_{y}<0$ region) of the
scattered FP intensity distribution as the intensity pattern for $k_{y}>0$
region is anti-symmetric with respect to $k_{y}<0$ region (Supplementary
information S5).

The detection mechanism is robust with respect to the inversion of the orbital
flow. We prove this by investigating the scattered FP intensity distribution
for $\sigma_{+}$ polarized and $\sigma_{-}$ polarized $\textrm{LG}^{-1}_{0}$
beam as shown in Figure 5(a) and (b) respectively. The corresponding intensity
distribution profile in sub critical region of the the $k_{y}>0$ half for the
$\sigma_{+}$ and $\sigma_{-}$ polarized case is shown in Figure 5 (c) and (d)
respectively. The measured directionality is found to be $D_{sub}=0.054$ for
$\sigma_{+}$ and $D_{sub}=-0.132$ for $\sigma_{-}$ polarized
$\textrm{LG}^{-1}_{0}$ beam. The preference of scattering in this region is
same as that of the $\sigma_{+}$ and $\sigma_{-}$ polarized
$\textrm{LG}^{1}_{0}$ beam, indicating the same SAM sign. Whereas, for the
super-critical region, the intensity distribution profile shown in Figure 5
(e) ($D_{super}=0.369$) and (f) ($D_{super}=0.114$), the preference is
opposite to that of the $\textrm{LG}^{1}_{0}$ case, indicating the sign
inversion of the topological charge and hence the orbital as well as total
transverse energy flow of the beam. Therefore, owing to the polarization and
transverse energy flow dependent interaction of the AgNW with the incident
beam, we are able to unambiguously detect OAM and SAM of the incident beam.
The accuracy of the simultaneous detection of the SAM and OAM state depends
upon the width of the nanowire used for the scattering process. We have
performed the experiments for nanowires having width $\approx 200$ nm, and the
results have been elaborated in supplementary information S6. In general, the
directionality in the sub- and super-critical region depends upon the incident
polarization state of the beam. For $\textrm{LG}^{1}_{0}$ beam, the elliptical
polarization dependence of the directionality and the reconstruction of
elliptical polarization states are investigated in supplementary information
S7.

Figure 6: Experimentally measured scattered FP intensity distribution of (a)
$\sigma_{+}$ polarized and (b) $\sigma_{-}$ polarized $\textrm{LG}^{2}_{0}$
beam. The corresponding intensity distribution profile in the sub-critical
region of the $k_{y}>0$ part of the FP for $\sigma_{+}$ and $\sigma_{-}$
polarized $\textrm{LG}^{2}_{0}$ beam is shown in (c) and (d) respectively. The
same for super-critical region is shown in (e) and (f). The black arrows
indicate the effective directionality.

###  2.5 Detection of OAM and SAM for higher order LG beam

Next, we explore the possibility of simultaneous detection of higher order OAM
as well as SAM and investigate the FP intensity distribution for higher order
LG beam: $\textrm{LG}^{2}_{0}$ beam. Figure 6(a) and (b) shows the FP
intensity distribution for $\sigma_{+}$ and $\sigma_{-}$ polarized
$\textrm{LG}^{2}_{0}$ beam respectively. Here we observe 3 lobes, instead of 2
like in the case of $\textrm{LG}^{1}_{0}$ beam, indicating the magnitude of
the topological charge ($n+1$ lobes for $\textrm{LG}^{n}_{0}$ beam)42. The
corresponding intensity distribution in the sub-critical region of the
$k_{y}>0$ half is shown in Figure 6 (c) and (d). The same for the super-
critical region is shown in figure 6(e) and (f). From both figure 6(c) and (d)
it can be seen that the central lobe remains of almost same intensity for
$\sigma_{+}$ and $\sigma_{-}$ polarized case, whereas the gray-shaded right-
most (numbered $1$) and left-most (numbered $3$) lobes undergo intensity
change. Thus we consider the integrated intensity of these two lobes only (1:
NA $0.5$ to $0.25$ and 2: NA $-0.25$ to $-0.5$). The directionality values for
the sub-critical region intensity distribution in figure 6 (c) and (d) are
$D_{sub}=0.029$ and $D_{sub}=-0.157$ respectively. As before, the sign
inversion of directionality allows us to accurately detect the SAM sign and
the signs are same as in the case of $\textrm{LG}^{1}_{0}$ beam. On the other
hand, the directionality for the super-critical region, $D_{super}=-0.278$ for
$\sigma_{+}$ polarized and $D_{super}=-0.555$ for $\sigma_{-}$ polarized case
indicates the OAM sign.

It is important to notice that the directionality values we get in the sub-
critical region for $\sigma_{+}$ polarized beam is lower for
$\textrm{LG}^{2}_{0}$ beam than that of $\textrm{LG}^{1}_{0}$ beam. This is
because, as the topological charge magnitude of LG beam increases, the
magnitude of the orbital flow increases (supplementary information S8), and
dominates the total energy flow and hence it will be more difficult to
distinguish the effect of the spin flow. Additionally, presence of $n+1$ lobes
for large $n$ for $\textrm{LG}^{n}_{0}$ beam scattering 42 in a small
scattered width ($k_{x}=$ $-0.5$ to $0.5$) in the FP intensity distribution,
will make it harder for us to detect the fluctuation in the polarization
dependent magnitude of the lobes.

##  3 Conclusion

To summarize, we have experimentally shown how forward scattering from a
single AgNW can be used for simultaneous unambiguous sensing of spin and
orbital angular momentum by employing Fourier plane microscopy. By analyzing
the intensity lobes due to circularly polarized LG beam in the sub- and super-
ciritical region in the Fourier plane, we show that the preference of
scattering in the sub-critical region changes when the spin is inverted,
whereas the super-critical region intensity distribution preference depends
predominantly on the orbital energy flow of the incident beam, thus allowing
us unambiguous detection of both. The incident linear polarization dependent
behaviour of the AgNW of scattering within the sub-critical region or both
sub- and super-critical region in the FP, allows the detection mechanism. The
focal optical field calculation gives us more insight onto the scattering
mechanism as well as the resulting directionality. Thus, our work emphasizes
spin-orbit interaction based simple single element chip-scale detection
configuration with prominent spin and orbital angular momentum detection
capability. Additionally, the spin-orbit interaction dependant scattering may
lead to opportunities in harnessing the resulting optical forces for optical
manipulation at nano scale as well as furthering our understanding of spin-
orbit interaction.

## Acknowledgements

This work was partially funded by Air Force Research Laboratory grant
(FA2386-18-1-4118 R&D 18IOA118), DST Energy Science grant (SR/NM/TP-13/2016)
and Swarnajayanti fellowship grant (DST/SJF/PSA-02/2017-18). DP and DKS
acknowledges Sunny Tiwari, Rahul Chand, Vandana Sharma, Shailendra K. Chaubey
and Chetna Taneja for fruitful discussions.

## References

  * Allen et al. 1992 Allen, L.; Beijersbergen, M. W.; Spreeuw, R. J. C.; Woerdman, J. P. Orbital angular momentum of light and the transformation of Laguerre-Gaussian laser modes. _Phys. Rev. A_ 1992, _45_ , 8185–8189
  * Molina-Terriza et al. 2007 Molina-Terriza, G.; Torres, J. P.; Torner, L. Twisted photons. _Nature Physics_ 2007, _3_ , 305–310
  * Franke-Arnold et al. 2008 Franke-Arnold, S.; Allen, L.; Padgett, M. Advances in optical angular momentum. _Laser & Photonics Reviews_ 2008, _2_ , 299–313
  * Yao and Padgett 2011 Yao, A. M.; Padgett, M. J. Orbital angular momentum: origins, behavior and applications. _Adv. Opt. Photon._ 2011, _3_ , 161–204
  * Shen et al. 2019 Shen, Y.; Wang, X.; Xie, Z.; Min, C.; Fu, X.; Liu, Q.; Gong, M.; Yuan, X. Optical vortices 30 years on: OAM manipulation from topological charge to multiple singularities. _Light: Science & Applications_ 2019, _8_ , 90
  * Bekshaev et al. 2011 Bekshaev, A.; Bliokh, K. Y.; Soskin, M. Internal flows and energy circulation in light beams. _Journal of Optics_ 2011, _13_ , 053001
  * Bliokh et al. 2010 Bliokh, K. Y.; Alonso, M. A.; Ostrovskaya, E. A.; Aiello, A. Angular momenta and spin-orbit interaction of nonparaxial light in free space. _Phys. Rev. A_ 2010, _82_ , 063825
  * Cardano and Marrucci 2015 Cardano, F.; Marrucci, L. Spin–orbit photonics. _Nature Photonics_ 2015, _9_ , 776–778
  * Aiello et al. 2015 Aiello, A.; Banzer, P.; Neugebauer, M.; Leuchs, G. From transverse angular momentum to photonic wheels. _Nature Photonics_ 2015, _9_ , 789–795
  * Bliokh et al. 2015 Bliokh, K. Y.; Rodríguez-Fortuño, F. J.; Nori, F.; Zayats, A. V. Spin–orbit interactions of light. _Nature Photonics_ 2015, _9_ , 796–808
  * Neugebauer et al. 2019 Neugebauer, M.; Nechayev, S.; Vorndran, M.; Leuchs, G.; Banzer, P. Weak Measurement Enhanced Spin Hall Effect of Light for Particle Displacement Sensing. _Nano Letters_ 2019, _19_ , 422–425
  * Zhao et al. 2007 Zhao, Y.; Edgar, J. S.; Jeffries, G. D. M.; McGloin, D.; Chiu, D. T. Spin-to-Orbital Angular Momentum Conversion in a Strongly Focused Optical Beam. _Phys. Rev. Lett._ 2007, _99_ , 073901
  * Nechayev et al. 2019 Nechayev, S.; Eismann, J. S.; Leuchs, G.; Banzer, P. Orbital-to-spin angular momentum conversion employing local helicity. _Phys. Rev. B_ 2019, _99_ , 075155
  * Zhao et al. 2009 Zhao, Y.; Shapiro, D.; Mcgloin, D.; Chiu, D. T.; Marchesini, S. Direct observation of the transfer of orbital angular momentum to metal particles from a focused circularly polarized Gaussian beam. _Opt. Express_ 2009, _17_ , 23316–23322
  * Rodríguez-Herrera et al. 2010 Rodríguez-Herrera, O. G.; Lara, D.; Bliokh, K. Y.; Ostrovskaya, E. A.; Dainty, C. Optical Nanoprobing via Spin-Orbit Interaction of Light. _Phys. Rev. Lett._ 2010, _104_ , 253601
  * Wang et al. 2012 Wang, J.; Yang, J.-Y.; Fazal, I. M.; Ahmed, N.; Yan, Y.; Huang, H.; Ren, Y.; Yue, Y.; Dolinar, S.; Tur, M.; Willner, A. E. Terabit free-space data transmission employing orbital angular momentum multiplexing. _Nature Photonics_ 2012, _6_ , 488–496
  * Willner et al. 2015 Willner, A. E. et al. Optical communications using orbital angular momentum beams. _Adv. Opt. Photon._ 2015, _7_ , 66–106
  * Lei et al. 2015 Lei, T.; Zhang, M.; Li, Y.; Jia, P.; Liu, G. N.; Xu, X.; Li, Z.; Min, C.; Lin, J.; Yu, C.; Niu, H.; Yuan, X. Massive individual orbital angular momentum channels for multiplexing enabled by Dammann gratings. _Light: Science & Applications_ 2015, _4_ , e257–e257
  * Xie et al. 2018 Xie, Z.; Lei, T.; Li, F.; Qiu, H.; Zhang, Z.; Wang, H.; Min, C.; Du, L.; Li, Z.; Yuan, X. Ultra-broadband on-chip twisted light emitter for optical communications. _Light: Science & Applications_ 2018, _7_ , 18001–18001
  * Mair et al. 2001 Mair, A.; Vaziri, A.; Weihs, G.; Zeilinger, A. Entanglement of the orbital angular momentum states of photons. _Nature_ 2001, _412_ , 313–316
  * Padgett et al. 1996 Padgett, M.; Arlt, J.; Simpson, N.; Allen, L. An experiment to observe the intensity and phase structure of Laguerre–Gaussian laser modes. _American Journal of Physics_ 1996, _64_ , 77–82
  * Leach et al. 2002 Leach, J.; Padgett, M. J.; Barnett, S. M.; Franke-Arnold, S.; Courtial, J. Measuring the Orbital Angular Momentum of a Single Photon. _Phys. Rev. Lett._ 2002, _88_ , 257901
  * Leach et al. 2004 Leach, J.; Courtial, J.; Skeldon, K.; Barnett, S. M.; Franke-Arnold, S.; Padgett, M. J. Interferometric Methods to Measure Orbital and Spin, or the Total Angular Momentum of a Single Photon. _Phys. Rev. Lett._ 2004, _92_ , 013601
  * Felde 2004 Felde, C. V. Diffraction diagnostics of phase singularities at optical beams. Sixth International Conference on Correlation Optics. 2004; pp 67 – 76
  * Jack et al. 2008 Jack, B.; Padgett, M. J.; Franke-Arnold, S. Angular diffraction. _New Journal of Physics_ 2008, _10_ , 103013
  * de Araujo and Anderson 2011 de Araujo, L. E. E.; Anderson, M. E. Measuring vortex charge with a triangular aperture. _Opt. Lett._ 2011, _36_ , 787–789
  * Bekshaev et al. 2020 Bekshaev, A.; Mikhaylovskaya, L.; Patil, S.; Kumar, V.; Singh, R. P. Optical-vortex diagnostics via Fraunhofer slit diffraction with controllable wavefront curvature. _J. Opt. Soc. Am. A_ 2020, _37_ , 780–786
  * Pu et al. 2015 Pu, M. et al. Catenary optics for achromatic generation of perfect optical angular momentum. _Science Advances_ 2015, _1_ , e1500396
  * Guo et al. 2016 Guo, Y.; Pu, M.; Zhao, Z.; Wang, Y.; Jin, J.; Gao, P.; Li, X.; Ma, X.; Luo, X. Merging Geometric Phase and Plasmon Retardation Phase in Continuously Shaped Metasurfaces for Arbitrary Orbital Angular Momentum Generation. _ACS Photonics_ 2016, _3_ , 2022–2029
  * Genevet et al. 2012 Genevet, P.; Lin, J.; Kats, M. A.; Capasso, F. Holographic detection of the orbital angular momentum of light with plasmonic photodiodes. _Nature Communications_ 2012, _3_ , 1278
  * Chen et al. 2018 Chen, J.; Chen, X.; Li, T.; Zhu, S. On-Chip Detection of Orbital Angular Momentum Beam by Plasmonic Nanogratings. _Laser & Photonics Reviews_ 2018, _12_ , 1700331
  * Zhang et al. 2020 Zhang, S.; Huo, P.; Zhu, W.; Zhang, C.; Chen, P.; Liu, M.; Chen, L.; Lezec, H. J.; Agrawal, A.; Lu, Y.; Xu, T. Broadband Detection of Multiple Spin and Orbital Angular Momenta via Dielectric Metasurface. _Laser & Photonics Reviews_ 2020, _14_ , 2000062
  * Sun et al. 2002 Sun, Y.; Yin, Y.; Mayers, B. T.; Herricks, T.; Xia, Y. Uniform Silver Nanowires Synthesis by Reducing AgNO3 with Ethylene Glycol in the Presence of Seeds and Poly(Vinyl Pyrrolidone). _Chemistry of Materials_ 2002, _14_ , 4736–4745
  * Visaveliya and Köhler 2013 Visaveliya, N.; Köhler, J. M. A self-seeding synthesis of Ag microrods of tuned aspect ratio: ascorbic acid plays a key role. _Nanotechnology_ 2013, _24_ , 345604
  * Yang et al. 2016 Yang, H.; Qiu, M.; Li, Q. Identification and control of multiple leaky plasmon modes in silver nanowires. _Laser & Photonics Reviews_ 2016, _10_ , 278–286
  * Taneja et al. 2021 Taneja, C.; Paul, D.; Pavan Kumar, G. V. Publisher’s Note: “Experimental observation of transverse spin of plasmon polaritons in a single crystalline silver nanowire” [Appl. Phys. Lett. 119, 161108 (2021)]. _Applied Physics Letters_ 2021, _119_ , 219901
  * Vasista et al. 2018 Vasista, A. B.; Jog, H.; Heilpern, T.; Sykes, M. E.; Tiwari, S.; Sharma, D. K.; Chaubey, S. K.; Wiederrecht, G. P.; Gray, S. K.; Kumar, G. V. P. Differential Wavevector Distribution of Surface-Enhanced Raman Scattering and Fluorescence in a Film-Coupled Plasmonic Nanowire Cavity. _Nano Letters_ 2018, _18_ , 650–655, PMID: 29244518
  * Tiwari et al. 2020 Tiwari, S.; Taneja, C.; Sharma, V.; Vasista, A. B.; Paul, D.; Pavan Kumar, G. V. Dielectric Microsphere Coupled to a Plasmonic Nanowire: A Self-Assembled Hybrid Optical Antenna. _Advanced Optical Materials_ 2020, _8_ , 1901672
  * Wei et al. 2011 Wei, H.; Wang, Z.; Tian, X.; Käll, M.; Xu, H. Cascaded logic gates in nanophotonic plasmon networks. _Nature Communications_ 2011, _2_ , 387
  * Lopez-Diaz et al. 2015 Lopez-Diaz, D.; Merino, C.; Velázquez, M. M. Modulating the Optoelectronic Properties of Silver Nanowires Films: Effect of Capping Agent and Deposition Technique. _Materials_ 2015, _8_ , 7622–7633
  * Sharma et al. 2018 Sharma, D. K.; Kumar, V.; Vasista, A. B.; Chaubey, S. K.; Kumar, G. V. P. Spin-Hall effect in the scattering of structured light from plasmonic nanowire. _Opt. Lett._ 2018, _43_ , 2474–2477
  * Sharma et al. 2019 Sharma, D. K.; Kumar, V.; Vasista, A. B.; Paul, D.; Chaubey, S. K.; Kumar, G. V. P. Optical Orbital Angular Momentum Read-Out Using a Self-Assembled Plasmonic Nanowire. _ACS Photonics_ 2019, _6_ , 148–153
  * Paul et al. 2021 Paul, D.; Sharma, D. K.; Kumar, G. V. P. Focused linearly-polarized-light scattering from a silver nanowire: Experimental characterization of the optical spin-Hall effect. _Phys. Rev. A_ 2021, _103_ , 013520
  * Wolf and Gabor 1959 Wolf, E.; Gabor, D. Electromagnetic diffraction in optical systems - I. An integral representation of the image field. _Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences_ 1959, _253_ , 349–357
  * Richards et al. 1959 Richards, B.; Wolf, E.; Gabor, D. Electromagnetic diffraction in optical systems, II. Structure of the image field in an aplanatic system. _Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences_ 1959, _253_ , 358–379
  * Novotny and Hecht 2006 Novotny, L.; Hecht, B. _Principles of Nano-Optics_ ; Cambridge University Press, 2006
  * Kerber et al. 2017 Kerber, R. M.; Fitzgerald, J. M.; Reiter, D. E.; Oh, S. S.; Hess, O. Reading the Orbital Angular Momentum of Light Using Plasmonic Nanoantennas. _ACS Photonics_ 2017, _4_ , 891–896
  * Kurvits et al. 2015 Kurvits, J. A.; Jiang, M.; Zia, R. Comparative analysis of imaging configurations and objectives for Fourier microscopy. _J. Opt. Soc. Am. A_ 2015, _32_ , 2082–2092
  * Vasista et al. 2019 Vasista, A. B.; Sharma, D. K.; Kumar, G. P. _digital Encyclopedia of Applied Physics_ ; American Cancer Society, 2019; pp 1–14

See SI_LPR_Paul


# Deterministic Epidemic Models For Ebola Infection With Time-dependent
Controls.

E. Okyere Department of Mathematics and Statistics, UENR, Sunyani,
<EMAIL_ADDRESS>J. D. Ankamah Department of Mathematics and
Statistics, UENR<EMAIL_ADDRESS>A. K. Hunkpe Department
of Mathematics and Statistics, UENR, Sunyani<EMAIL_ADDRESS>D.
Mensah Department of Mathematics and Statistics, UENR, Sunyani,
<EMAIL_ADDRESS>
###### Abstract

In this paper, we have studied epidemiological models for Ebola infection
using nonlinear ordinary differential equations and optimal control theory. We
considered optimal control analysis of SIR and SEIR models for the deadly
Ebola infection using vaccination, treatment and educational campaign as time-
dependent controls functions. We have applied indirect methods to study
existing deterministic optimal control epidemic models for Ebola virus
disease. These methods in optimal control are based on Hamiltonian function
and the Pontryagin’s maximum principle to construct adjoint equations and
optimality systems. The forward-backward sweep numerical scheme with fourth-
order Runge-Kutta method is used to solve the optimality system for the
various control strategies. From our numerical illustrations, we can conclude
that, effective educational campaigns and vaccination of susceptible
individuals as were as effective treatments of infected individuals can help
reduce the disease transmission.

##  1 Introduction

The re-emergence of the Ebola virus disease in $2014-2016$ in West Africa has
been classified as the largest outbreak since the disease was first discovered
in DRC in 1976 (WHO, 2019). This highly infectious and deadly disease has
claimed many lives and caused huge economic burden in the affected West Africa
Countries. The recent $2018-2019$ outbreak in eastern DRC is a very complex
situation due to insecurities which is seriously affecting public health
workers response activities (WHO, 2019).

Mathematical modelling of epidemics has contributed significantly in
understanding the dynamical behaviour and control of infectious diseases
(Hethcote, 2000). The complex dynamics of Ebola virus disease has attracted
the attention of many researchers who are interested in epidemiological
modeling [see, e.g, Tulu et al. (2017), Dong et al. (2015), Althaus (2014),
Chowell and Nishiura (2014), Ngwa and Teboh-Ewungkem (2016), Weitz and Dushoff
(2015), Ndanguza et al. (2017), Diaz et al. (2018), Atangana and Goufo (2014),
Jiang et al. (2017), Berge et al. (2017), Xia et al. (2015), Khan et al.
(2015), Chowell et al. (2004), Goufo et al. (2016), Agusto et al. (2015), Luo
et al. (2019), Dénes and Gumel (2019)].

System of nonlinear equations that incorporates optimal control dynamics are
key mathematical tools that are used in compartmental modelling to understand
the spread of infectious diseases. In Sharomi and Malik (2017), they conducted
a comprehensive survey on optimal control modeling of several infectious
diseases. Zakary et al. (2017) formulated epidemic model for controlling Ebola
outbreak. A mathematical model for Ebola disease with time-dependent controls
is proposed and numerically analyzed in (Bonyah et al., 2016). The authors in
(Ahmad et al., 2016) proposed and studied optimal control epidemic model for
Ebola virus infection. Area et al. (2018) introduced and studied deterministic
epidemic model for Ebola disease that incorporates optimal control dynamics
with vaccination control. The authors in (Grigorieva et al., 2017) developed
and analyzed controlled dynamical model for Ebola virus infection.

Our present work is motivated by indirect methods in optimal control theory.
These methods are based on the Pontryagin’s maximum principle and Hamiltonian
function to construct adjoint equations and optimality systems for optimal
control problems. Indirect methods have widely been applied by several
authors, see the detailed survey on optimal modeling of infectious diseases by
the authors in (Sharomi and Malik, 2017) and references therein. Recently,
indirect methods in optimal control have also been used to study the dynamical
behaviour and control of Zika virus disease [see, e.g, Bonyah et al. (2017),
Khan et al. (2019), Miyaoka et al. (2019), Bonyah et al. (2019)]. We are also
inspired by the forward-backward numerical scheme with fourth order Runge-
Kutta method for optimal control problems described in (Lenhart and Workman,
2007).

In the works by the authors in (Rachah and Torres, 2015, 2016), they
formulated and analyzed SIR and SEIR optimal control models for Ebola
infection using direct methods. In both papers, they formulated objective
functionals with dynamical state constrained equations. In (Rachah and Torres,
2015), they proposed one optimal control strategy for the SIR model and the
same authors in (Rachah and Torres, 2016) proposed three different optimal
control strategies for the SEIR model. After their models formulations, they
applied the ACADO solver which is an automatic control and dynamic
optimization tool to perform their numerical simulations. In our present work,
will apply indirect methods to study the mathematical models proposed in
(Rachah and Torres, 2015, 2016). In this study, we will also propose one
additional optimal control strategy for the SIR model. We will analytically
construct Hamiltonian function and optimality system for the various optimal
control strategies. We will then apply the forward-backward sweep method with
fourth-order Runge-Kutta method to perform numerical simulations in Matlab.

The rest of the paper is organized as follows. In section 2, we will introduce
the classic SIR epidemic model to describe the dynamical behaviour of Ebola
infection. We will then consider optimal control problems for the SIR model
with two control strategies. In section 3, we will present the basic SEIR
epidemic model to describe the transmission dynamics of Ebola virus disease.
We will further consider SEIR optimal control problems for this infectious
disease with three different control strategies. In all these sections, we
will formulate Hamiltonian functions and then apply the Pontryagin’s maximum
principle to construct adjoint equations and optimality system for the various
optimal control strategies. We will also perform numerical simulations for the
optimality systems. Finally, we will conclude the paper in section 4.

##  2 SIR Model

In this section, we introduce the classic SIR epidemic model to describe the
dynamical behaviour of Ebola infection. The model assumes constant population
size with no vital dynamics (Hethcote, 2000). The population is divided into
three different classes with $S(t),\ I(t)$ and $R(t)$ representing
Susceptible, Infected and Recovered individuals respectively. As in (Rachah
and Torres, 2015), the nonlinear dynamical system describing Ebola infection
is given by

$\begin{cases}\frac{dS(t)}{dt}={-\frac{\nu S(t)I(t)}{N}},\ \ \qquad\qquad
S(0)=S_{0}\geq 0,&\\\ \\\ \frac{dI(t)}{dt}={\frac{\nu
S(t)I(t)}{N}}-{\delta}{I(t)},\qquad I(0)=I_{0}\geq 0,\\\ \\\
\frac{dR(t)}{dt}={\delta I(t)},\ \qquad\qquad\qquad R(0)=R_{0}\geq
0.\end{cases}$ (1)

with $S(t)+I(t)+R(t)=N,$

where $\nu$ is the infection rate and $\delta$ is the recovery rate.

In this study, we will work with proportional quantities instead of the actual
populations by scaling each state variable ($S,\ I,\ R$) by the total
population. For this purpose, we introduce new state variables ($s,\ i,\ r$)
which are expressed in terms of the original state variables and the total
population given as follows;

$\begin{cases}s(t)=\frac{S(t)}{N},\\\ \\\ i(t)=\frac{I(t)}{N},\\\ \\\
r(t)=\frac{R(t)}{N}.\end{cases}$ (2)

Differentiating the scaling equation (2) with respect to time $t$ and using
the unscaled model problem (1), we obtain the scaled SIR model describing the
transmission dynamics of Ebola is as follows

$\begin{cases}\frac{ds(t)}{dt}={-\nu s(t)i(t)},\qquad\qquad s(0)=s_{0}\geq
0,&\\\ \\\ \frac{di(t)}{dt}={\nu s(t)i(t)}-{\delta i(t)},\qquad i(0)=i_{0}\geq
0,\\\ \\\ \frac{dr(t)}{dt}={\delta i(t)},\ \ \qquad\qquad\qquad r(0)=r_{0}\geq
0.\end{cases}$ (3)

with $s(t)+i(t)+r(t)=1$,

where the new state variables $s(t),\ i(t)$ and $r(t)$ represent proportions
of Susceptible, Infected and Recovered individuals respectively.

As in (Rachah and Torres, 2015), we will use this scaled SIR model (3) in
subsections 2.1 and 2.2 to formulate optimal control models for Ebola
infection with two different optimal control strategies.

###  2.1 SIR Model with Optimal Control Strategy 1

In this subsection, we consider a controlled dynamical system for the scaled
SIR model (3), using vaccination as time dependent control function,
$\eta(t)$. For this optimal control problem, we want to reduce the number of
infected individuals and the cost of vaccination. Therefore, in this control
strategy, we minimize the objective functional $J(\eta)$ given by

$J(\eta)=\int^{t_{f}}_{0}\left[i(t)+\dfrac{B}{2}\eta^{2}(t)\right]dt$ (4)

subject to:

$\begin{cases}\frac{ds(t)}{dt}={-\nu s(t)i(t)-\eta(t)s(t)},\ \qquad
s(0)=s_{0}\geq 0,\\\ \\\ \frac{di(t)}{dt}={\nu s(t)s(t)-\delta i(t)},\
\qquad\qquad i(0)=i_{0}\geq 0,\\\ \\\ \frac{dr(t)}{dt}={\delta
i(t)+\eta(t)s(t)},\ \ \qquad\qquad r(0)=r_{0}\geq 0.\end{cases}$ (5)

where the control set is given as:

${\mathcal{G}}_{1}=\\{\eta\mathrel{\mathop{\mathchar 58\relax}}\eta(t)\
\textrm{is lebesgue measurable},\ 0\leq\eta(t)\leq 1,\ t\in[0,t_{f}]\\}$

and the positive parameter $A$ is the weight on cost of vaccination.

The Hamiltonian function corresponding to the objective functional (4) and the
dynamic constrained state system (5) is given by

$\begin{split}H&=i(t)+\dfrac{B}{2}\eta^{2}(t)\\\ \\\ &+\varphi_{1}\left[-\nu
s(t)i(t)-\eta(t)s(t)\right]\\\ \\\ &+\varphi_{2}\left[\nu s(t)i(t)-\delta
i(t)\right]\\\ \\\ &+\varphi_{3}\left[\delta
i(t)+\eta(t)s(t)\right]\end{split}$ (6)

We then apply the Maximum Principle proposed by authors in (Pontryagin et al.,
1962) to determine an optimal solution as follows:

Suppose that $(w,\eta)$ is an optimal solution for a controlled dynamical
system, then there exist adjoint vector function $\varphi=(\varphi_{1},\
\varphi_{2},\ \cdots\varphi_{n})$ which satisfy the system below;

$\begin{cases}\frac{dw}{dt}={\frac{\partial
H(t,w,\eta,\varphi)}{\partial\varphi}},\\\ \\\ 0={\frac{\partial
H(t,w,\eta,\varphi)}{\partial\eta}},\\\ \\\
\frac{d\varphi}{dt}=-{\frac{\partial H(t,w,\eta,\varphi)}{\partial
w}}.\end{cases}$ (7)

Following the constructed Hamiltonian function (6) and equation (7), we
present the adjoint equations and the control characterisation as follows

###### Theorem 1

Let ${\eta}^{*}$ be an optimal control and optimal state solutions $s^{*},\
i^{*},\ r^{*}$ of the corresponding controlled dynamical system (4)-(5) that
minimize $J(\eta)$ over ${\mathcal{G}}_{1}$. Then there exist adjoint
variables $\varphi_{i}$ for $i=1,\ 2,3$ which satisfy the system below

$\begin{cases}\displaystyle\frac{d\varphi_{1}}{dt}=\displaystyle{\nu
i^{*}(t)(\varphi_{1}-\varphi_{2})+\eta^{*}(t)(\varphi_{1}-\varphi_{3})},\qquad\\\
\\\ \displaystyle\frac{d\varphi_{2}}{dt}=\displaystyle{-1+\nu
s^{*}(t)(\varphi_{1}-\varphi_{2})+\delta(\varphi_{2}-\varphi_{3})},\qquad\\\
\\\ \displaystyle\frac{d\varphi_{3}}{dt}=\displaystyle{0}.\end{cases}$ (8)

with transversality conditions $\varphi_{1}(t_{f})=0,\ \varphi_{2}(t_{f})=0,\
\varphi_{3}(t_{f})=0$

and the control ${\eta}^{*}$ satisfies the optimality condition.

$\displaystyle{\eta}^{*}(t)=min\left\\{max\left\\{0,\dfrac{s^{*}(t)}{B}(\varphi_{1}-\varphi_{3})\right\\},1\right\\}$
(9)

Proof. To derive the system of equations describing the adjoint variables and
transversality conditions, we apply the Pontryagin’s Maximum Principle and the
Hamiltonian function (6) as follows

$\begin{cases}\frac{d\varphi_{1}}{dt}=-\frac{\partial H}{\partial s},\\\ \\\
\frac{d\varphi_{2}}{dt}=-\frac{\partial H}{\partial i},\\\ \\\
\frac{d\varphi_{3}}{dt}=-\frac{\partial H}{\partial r}.\end{cases}$ (10)

with

$\varphi_{i}(t_{f})=0,\ \ i=1,\ 2,\ 3.$ (11)

By using the property of the control space ${\mathcal{G}}_{1}$ and solving the
differential equation given by

$\frac{\partial H}{\partial\eta}=0$ (12)

on the interior of the control set, the characterization of the optimal
control is obtained as given by equation (9).

###  2.2 SIR Model with Optimal Control Strategy 2

In this subsection, we present a controlled dynamical system by incorporating
treatment control, $\eta_{1}(t)$ and educational campaign control,
$\eta_{2}(t)$ into the scaled SIR model (3). For this optimal control problem,
we want to reduce the number of infected individuals, the cost of educational
campaign and cost of treatment. Therefore, in this control strategy, we
minimise the objective functional given as

$J(\eta_{1},\eta_{2})=\int^{t_{f}}_{0}\left[C_{1}i(t)+\dfrac{C_{2}}{2}\eta^{2}_{1}(t)+\dfrac{C_{3}}{2}\eta^{2}_{2}(t)\right]dt$
(13)

subject to

$\begin{cases}\frac{ds(t)}{dt}={-\nu
s(t)i(t)-\eta_{2}(t)s(t)},\qquad\qquad\qquad s(0)=s_{0}\geq 0,&\\\ \\\
\frac{di(t)}{dt}={\nu s(t)i(t)}-{\delta i(t)-\eta_{1}(t)i(t)},\qquad\qquad
i(0)=i_{0}\geq 0,\\\ \\\ \frac{dr(t)}{dt}={\delta
i(t)+\eta_{1}(t)i(t)+\eta_{2}(t)s(t)}\qquad\qquad r(0)=r_{0}\geq
0.\end{cases}$ (14)

where the control set is given by

${\mathcal{G}}_{2}=\\{(\eta_{1},\eta_{2})\mathrel{\mathop{\mathchar
58\relax}}\eta_{i}(t)\ \textrm{is lebesgue measurable},\ 0\leq\eta_{i}(t)\leq
1,\ t\in[0,t_{f}],\ \textrm{for}\ i=1,2\\}$

and $C_{1}$ is a positive constant to keep a balance in the size of $i(t)$.
The positive parameters $C_{2}$ and $C_{3}$ are the weight on cost for
treatments and educational campaigns.

The Hamiltonian function corresponding to equations (13) and (14) is given as

$\begin{split}H&=C_{1}i(t)+\dfrac{C_{2}}{2}\eta^{2}_{1}(t)+\dfrac{C_{3}}{2}\eta^{2}_{2}(t)\\\
\\\ &+\varphi_{1}\left[-\nu s(t)i(t)-\eta_{2}(t)s(t)\right]\\\ \\\
&+\varphi_{2}\left[\nu s(t)i(t)-\delta i(t)-\eta_{1}(t)i(t)\right]\\\ \\\
&+\varphi_{3}\left[\delta
i(t)+\eta_{1}(t)i(t)+\eta_{2}(t)s(t)\right]\end{split}$ (15)

We then apply the Maximum Principle proposed by authors in (Pontryagin et al.,
1962) to determine an optimal solution as

Assume that $(w,\eta)$ is an optimal solution for a controlled dynamical
system, then there exist adjoint vector function $\varphi=(\varphi_{1},\
\varphi_{2},\ \cdots\varphi_{n})$ which satisfy the following equations

$\begin{cases}\frac{dw}{dt}={\frac{\partial
H(t,w,\eta,\varphi)}{\partial\varphi}},\\\ \\\ 0={\frac{\partial
H(t,w,\eta,\varphi)}{\partial\eta}},\\\ \\\
\frac{d\varphi}{dt}=-{\frac{\partial H(t,w,\eta,\varphi)}{\partial
w}}.\end{cases}$ (16)

Using the formulated Hamiltonian function (6) and equation (16), the adjoint
equations and the control characterisation are presented in the following
Theorem.

###### Theorem 2

Let $\eta^{*}_{1}$ and $\eta^{*}_{2}$ be optimal control pair and optimal
state solutions $s^{*},\ i^{*},\ r^{*}$ of the corresponding controlled
dynamical system (13)-(14) that minimize $J(\eta_{1},\eta_{2})$ over
${\mathcal{G}}_{2}$. Then there exist adjoint variables $\varphi_{i}$ for
$i=1,\ 2,3$ satisfying

$\begin{cases}\displaystyle\frac{d\varphi_{1}}{dt}=\displaystyle{\nu
i^{*}(t)(\varphi_{1}-\varphi_{2})+\eta^{*}_{2}(t)(\varphi_{1}-\varphi_{3})},\qquad\\\
\\\ \displaystyle\frac{d\varphi_{2}}{dt}=\displaystyle{-C_{1}+\nu
s^{*}(t)(\varphi_{1}-\varphi_{2})+\delta(\varphi_{2}-\varphi_{3})+\eta^{*}_{1}(t)(\varphi_{2}-\varphi_{3})},\qquad\\\
\\\ \displaystyle\frac{d\varphi_{3}}{dt}=\displaystyle{0}.\end{cases}$ (17)

with transversality conditions $\varphi_{1}(t_{f})=0,\ \varphi_{2}(t_{f})=0,\
\varphi_{3}(t_{f})=0$

and the control pair $\eta^{*}_{1}$ and $\eta^{*}_{2}$ satisfies the
optimality conditions

$\begin{cases}\displaystyle\eta^{*}_{1}(t)=min\left\\{max\left\\{0,\dfrac{i^{*}(t)}{C_{2}}(\varphi_{2}-\varphi_{3})\right\\},1\right\\}\\\
\\\
\displaystyle\eta^{*}_{2}(t)=min\left\\{max\left\\{0,\dfrac{s^{*}(t)}{C_{3}}(\varphi_{1}-\varphi_{3})\right\\},1\right\\}\end{cases}$
(18)

###### Remark 1

We remark that, the proof of Theorem 2 and subsequent Theorems in the next
section of this study are similar to that of Theorem 1 and we therefore omit
their proofs.

###  2.3 Numerical Simulations and Discussions

This subsection deals with numerical solutions of optimal control problems
formulated in subsections 2.1 and 2.2. Analytical solution of nonlinear system
of differential equations with optimal control is a very hard task in
mathematical modelling. Therefore, in this study, we apply the efficient
forward-backward sweep numerical scheme with fourth-order Runge-Kutta method
to solve the optimality systems formulated for the two optimal control
strategies. This numerical scheme has extensively been described by the
authors in (Lenhart and Workman, 2007) in their mathematical modelling
textbook which is concern with the applications of differential equations and
optimal control theory. For our numerical simulations, we have adapted the
same model parameter values ($\nu=0.2,\ \delta=0.1$) as in the work by the
authors in (Rachah and Torres, 2015) and initial conditions: $s_{0}=0.95,\
i_{0}=0.05,\ r_{0}=0$. For the positive parameters in the objectives
functionals, we have assumed that $D=1,\ C_{1}=1,\ C_{2}=5$ and $C_{3}=5.$ In
Fig. 1, it is clear that, there is a significant decrease in the number of
susceptible individuals with control strategies 1, 2 than without control.
Figure 2 shows solution paths for infected individuals with control strategies
1, 2 and without control. This plot shows significant decrease in the infected
individuals with control strategies 1, 2 than without control. In Fig 3, there
is a rapid increase in recovered individuals with control strategies 1, 2 than
without control. Figures 4 and 5 represents control functions for optimal
control strategies 1 and 2 respectively.

Figure 1: Solution paths for susceptible individuals with two control
strategies and without control, $\nu=0.2$, $\delta=0.1$ Figure 2: Solution
paths for Infected individuals with two control strategies and without
control, $\nu=0.2$, $\delta=0.1$ Figure 3: Solution paths for Recovered
individuals with two control strategies and without control, $\nu=0.2$,
$\delta=0.1$

Figure 4: Optimal control function $\eta(t)$ for strategy 1 Figure 5: Optimal
control functions $\eta_{1}(t)$ and $\eta_{2}(t)$ for strategy 2

##  3 SEIR Model

In this section, we introduce the basic SEIR mathematical model to describe
the transmission dynamics of Ebola infection. This model also assumes constant
population size with no vital dynamics(birth and death rates). The total
population is divided into four different classes with $S(t),\ E(t),\ E(t)$
and $R(t)$ representing Susceptible, Exposed, Infectious and Recovered
individuals respectively. Following the compartmental modeling concepts in
(Hethcote, 2000), the nonlinear dynamical system describing Ebola virus
disease is given by

$\begin{cases}\displaystyle\frac{dS(t)}{dt}=\displaystyle{-\frac{\nu
S(t)I(t)}{N}},\ \ \ \qquad\qquad\qquad S(0)=S_{0}\geq 0,\\\ \\\
\displaystyle\frac{dE(t)}{dt}=\displaystyle{\frac{\nu
S(t)I(t)}{N}}-{\rho}{E(t)},\ \ \ \ \ \qquad E(0)=E_{0}\geq 0,\\\ \\\
\displaystyle\frac{dI(t}{dt}=\displaystyle{\rho E(t)}-{\delta}{I(t)},\ \
\qquad\qquad\qquad I(0)=I_{0}\geq 0,\\\ \\\
\displaystyle\frac{dR(t)}{dt}=\displaystyle{\delta
I(t)},\qquad\qquad\qquad\qquad\qquad R(0)=R_{0}\geq 0.\end{cases}$ (19)

with $S(t)+E(t)+I(t)+R(t)=N,$

where $\nu$ is the transmission rate, $\rho$ represent the infectious rate and
$\delta$ is the recovery rate.

As we did in section 2, we will work with proportional quantities instead of
the actual populations by scaling each state variable ($S(t),\ E(t),\ I(t),\
R(t)$) by the total population. We therefore introduce new state variables
($s(t),\ e(t),\ i(t),\ r(t)$) which are expressed in terms of the original
state variables and the total population is given by

$\begin{cases}s(t)=\frac{S(t)}{N},\\\ \\\ e(t)=\frac{E(t)}{N},\\\ \\\
i(t)=\frac{I(t)}{N},\\\ \\\ r(t)=\frac{R(t)}{N}.\end{cases}$ (20)

By differentiating the scaling equation (20) with respect to time $t$ and
using the unscaled model problem (19), we obtain the scaled SEIR model as
follows

$\begin{cases}\displaystyle\frac{ds(t)}{dt}=\displaystyle{-\nu s(t)i(t)},\
\qquad\qquad\qquad s(0)=s_{0}\geq 0,\\\ \\\
\displaystyle\frac{de(t)}{dt}=\displaystyle{\nu
s(t)i(t)}-{\rho}{e(t)},\qquad\qquad e(0)=e_{0}\geq 0,\\\ \\\
\displaystyle\frac{di(t)}{dt}=\displaystyle{\rho
e(t)}-{\delta}{i(t)},\qquad\qquad\qquad i(0)=i_{0}\geq 0,\\\ \\\
\displaystyle\frac{dr(r)}{dt}=\displaystyle{\delta i(t)},\ \ \
\qquad\qquad\qquad\qquad r(0)=r_{0}\geq 0.\end{cases}$ (21)

with $s(t)+e(t)+i(t)+r(t)=1,$

where the new state variables $s(t)\,e(t),\ i(t)$ and $r(t)$ represent
proportions of Susceptible, Infected and Recovered individuals respectively.

As in (Rachah and Torres, 2016), we will use this scaled SEIR model problem
(21) in subsections 3.1, 3.2 and 3.3 to construct optimal control models for
Ebola infection with three different control strategies.

###  3.1 SEIR Model with Optimal Control Strategy 1

In this subsection, we consider a controlled dynamical system for the scaled
SEIR model (21), using vaccination as time dependent control function,
$\eta(t)$. As we did in subsection 2.1, our main objective for this optimal
control problem is that, we want to reduce the number of infected individuals
and the cost of vaccination. Therefore in this control strategy, we minimize
the objective functional $J(\eta)$ given by

$J(\eta)=\int^{t_{f}}_{0}\left[i(t)+\dfrac{D}{2}\eta^{2}(t)\right]dt$ (22)

subject to:

$\begin{cases}\displaystyle\frac{ds(t)}{dt}=\displaystyle{-\nu
s(t)i(t)-\eta(t)s(t)},\ \qquad\qquad s(0)=s_{0}\geq 0,\\\ \\\
\displaystyle\frac{de(t)}{dt}=\displaystyle{\nu
s(t)i(t)}-{\rho}{e(t)},\qquad\qquad\qquad e(0)=e_{0}\geq 0,\\\ \\\
\displaystyle\frac{di(t)}{dt}=\displaystyle{\rho
e(t)}-{\delta}{i(t)},\qquad\qquad\qquad\qquad i(0)=i_{0}\geq 0,\\\ \\\
\displaystyle\frac{dr(t)}{dt}=\displaystyle{\delta i(t)+\eta(t)s(t)},\ \
\qquad\qquad\qquad r(0)=r_{0}\geq 0.\end{cases}$ (23)

where the control set is given by

${\mathcal{G}}_{3}=\\{\eta\mathrel{\mathop{\mathchar 58\relax}}\eta(t)\
\textrm{is lebesgue measurable},\ 0\leq\eta(t)\leq 1,\ t\in[0,t_{f}]\\}$

and the positive constant $D$ is the weight on the cost of vaccination.

The Hamiltonian function $H$ is then given as:

$\begin{split}H&=i(t)+\dfrac{D}{2}\eta^{2}(t)\\\ \\\ &+\varphi_{1}\left[-\nu
s(t)i(t)-\eta(t)s(t)\right]\\\ \\\ &+\varphi_{2}\left[\nu s(t)i(t)-\rho
e(t)\right]\\\ \\\ &+\varphi_{3}\left[\rho e(t)-\delta i(t)\right]\\\ \\\
&+\varphi_{4}\left[\delta i(t)+\eta(t)s(t)\right]\end{split}$ (24)

We then apply the Maximum Principle proposed by authors in (Pontryagin et al.,
1962) to determine an optimal solution as

Given that $(w,\eta)$ is an optimal solution for a controlled dynamical
system, then there exist adjoint vector function $\varphi=(\varphi_{1},\
\varphi_{2},\ \cdots\varphi_{n})$ which satisfy the following equations

$\begin{cases}\frac{dw}{dt}={\frac{\partial
H(t,w,\eta,\varphi)}{\partial\varphi}},\\\ \\\ 0={\frac{\partial
H(t,w,\eta,\varphi)}{\partial\eta}},\\\ \\\
\frac{d\varphi}{dt}=-{\frac{\partial H(t,w,\eta,\varphi)}{\partial
w}}.\end{cases}$ (25)

Using the formulated Hamiltonian function (24) and equation (25), the adjoint
equations and the control characterisation are presented in the following
Theorem.

###### Theorem 3

Let ${\eta}^{*}$ be an optimal control and optimal state solutions
$s^{*},e^{*},i^{*},r^{*}$ of the corresponding controlled dynamical system
(22)-(23) that minimize $J(\eta)$ over ${\mathcal{G}}_{3}$. Then there exist
adjoint variables $\varphi_{i}$ for $i=1,2,3,4$ satisfying

$\begin{cases}\displaystyle\frac{d\varphi_{1}}{dt}=\displaystyle{\nu
i^{*}(t)(\varphi_{1}-\varphi{2})+\eta^{*}(t)(\varphi_{1}-\varphi_{4})},\qquad\\\
\\\
\displaystyle\frac{d\varphi_{2}}{dt}=\displaystyle{\rho(\varphi_{2}-\varphi_{3})},\qquad\\\
\\\ \displaystyle\frac{d\varphi_{3}}{dt}=\displaystyle{-1+\nu
s^{*}(t)(\varphi_{1}-\varphi_{2})+\delta(\varphi_{3}-\varphi_{4})},\qquad\\\
\\\ \displaystyle\frac{d\varphi_{4}}{dt}=\displaystyle{0}.\qquad\end{cases}$
(26)

with transversality conditions $\varphi_{1}(t_{f})=0,\ \varphi_{2}(t_{f})=0,\
\varphi_{3}(t_{f})=0,\ \varphi_{4}(t_{f})=0$ and the control $\eta^{*}(t)$
satisfies the optimality condition.

$\displaystyle\eta^{*}(t)=min\left\\{max\left\\{0,\dfrac{s^{*}}{D}(\varphi_{1}-\varphi_{4})\right\\},1\right\\}$
(27)

###  3.2 SEIR Model with Optimal Control Strategy 2

In this subsection, we present a controlled dynamical system for the scaled
SEIR model (21), using vaccination as time dependent control function,
$\eta(t)$. Our main objective for this strategy is to reduce number of
infected and exposed individuals as were as the cost of vaccination. Therefore
in this control strategy, we minimize the objective functional $J(\eta)$ given
by

$J(\eta)=\int^{t_{f}}_{0}\left[K_{1}e(t)+K_{2}i(t)+\dfrac{K_{3}}{2}\eta^{2}(t)\right]dt$
(28)

subject to:

$\begin{cases}\displaystyle\frac{ds(t)}{dt}=\displaystyle{-\nu
s(t)i(t)-\eta(t)s(t)},\ \qquad\qquad s(0)=s_{0}\geq 0,\\\ \\\
\displaystyle\frac{de(t)}{dt}=\displaystyle{\nu s(t)i(t)}-{\rho}{e(t)},\ \ \ \
\ \ \qquad\qquad e(0)=e_{0}\geq 0,\\\ \\\
\displaystyle\frac{di(t)}{dt}=\displaystyle{\rho
e(t)}-{\delta}{i(t)},\qquad\qquad\qquad\qquad i(0)=i_{0}\geq 0,\\\ \\\
\displaystyle\frac{dr(t)}{dt}=\displaystyle{\delta i(t)+\eta(t)s(t)},\ \
\qquad\qquad\qquad r(0)=r_{0}\geq 0.\end{cases}$ (29)

where the control set is given as:

${\mathcal{G}}_{4}=\\{\eta\mathrel{\mathop{\mathchar 58\relax}}\eta(t)\
\textrm{is lebesgue measurable},\ 0\leq\eta(t)\leq 1,\ t\in[0,t_{f}]\\}$

and the parameters $K_{1}$ and $K_{2}$ are positive constants associated with
exposed individuals and infected individuals and $K_{3}$ is a positive weight
constant for control function $\eta(t)$.

The Hamiltonian $H$ is given by

$\begin{split}H&=K_{1}e(t)+K_{2}i(t)+\dfrac{K_{3}}{2}\eta^{2}(t)\\\ \\\
&+\varphi_{1}\left[-\nu s(t)i(t)-\eta(t)s(t)\right]\\\ \\\
&+\varphi_{2}\left[\nu s(t)i(t)-\rho e(t)\right]\\\ \\\
&+\varphi_{3}\left[\rho e(t)-\delta i(t)\right]\\\ \\\
&+\varphi_{4}\left[\delta i+\eta(t)s(t)\right]\end{split}$ (30)

We then apply the Maximum Principle proposed by authors in (Pontryagin et al.,
1962) to determine an optimal solution as

Assume that $(w,\eta)$ is an optimal solution for a controlled dynamical
system, then there exist adjoint vector function $\varphi=(\varphi_{1},\
\varphi_{2},\ \cdots\varphi_{n})$ which satisfy the following equations

$\begin{cases}\frac{dw}{dt}={\frac{\partial
H(t,w,\eta,\varphi)}{\partial\varphi}},\\\ \\\ 0={\frac{\partial
H(t,x,\eta,\varphi)}{\partial\eta}},\\\ \\\
\frac{d\varphi}{dt}=-{\frac{\partial H(t,w,\eta,\varphi)}{\partial
w}}.\end{cases}$ (31)

Using the formulated Hamiltonian function (30) and equation (31), the adjoint
equations and the control characterisation are given in the following Theorem
as

###### Theorem 4

Let ${\eta}^{*}$ be an optimal control and optimal state solutions
$s^{*},e^{*},i^{*},r^{*}$ of the corresponding controlled dynamical system
(28)-(29) that minimize $J(\eta)$ over ${\mathcal{G}}_{3}$. Then there exist
adjoint variables $\varphi_{i}$ for $i=1,2,3,4$ satisfying

$\begin{cases}\displaystyle\frac{d\varphi_{1}}{dt}=\displaystyle{\nu
i^{*}(t)(\varphi_{1}-\varphi{2})+\eta^{*}(t)(\varphi_{1}-\varphi_{4})},\qquad\\\
\\\
\displaystyle\frac{d\varphi_{2}}{dt}=\displaystyle{-K_{1}+\rho(\varphi_{2}-\varphi_{3})},\qquad\\\
\\\ \displaystyle\frac{d\varphi_{3}}{dt}=\displaystyle{-K_{2}+\nu
s^{*}(t)(\varphi_{1}-\varphi_{2})+\delta(\varphi_{3}-\varphi_{4})},\qquad\\\
\\\ \displaystyle\frac{d\varphi_{4}}{dt}=\displaystyle{0}.\end{cases}$ (32)

with transversality conditions $\varphi_{1}(t_{f})=0,\ \varphi_{2}(t_{f})=0,\
\varphi_{3}(t_{f})=0,\ \varphi_{4}(t_{f})=0$ and the control function
$\eta^{*}$ is given by

$\displaystyle\eta^{*}(t)=min\left\\{max\left\\{0,\dfrac{s^{*}(t)}{K_{3}}(\varphi_{1}-\varphi_{4})\right\\},1\right\\}$
(33)

###  3.3 SEIR Model with Optimal Control Strategy 3

In this subsection, we present a controlled dynamical problem by incorporating
treatment control, $\eta_{1}(t)$ and educational campaign control,
$\eta_{2}(t)$ into the SEIR model (3). As we presented in subsection 2.2, our
main aim for this controlled problem is to minimize the number of infected
individuals, the cost of educational campaigns and treatment. Therefore in
this control strategy, we minimise the objective functional given by

$J(\eta_{1},\eta_{2})=\int^{t_{f}}_{0}\left[D_{1}i(t)]+\dfrac{D_{2}}{2}\eta^{2}_{1}(t)+\dfrac{D_{3}}{2}\eta^{2}_{2}(t)\right]dt$
(34)

Subject to:

$\begin{cases}\displaystyle\frac{ds(t)}{dt}=\displaystyle{-\nu
s(t)i(t)-\eta_{2}(t)s(t)},\ \ \ \ \qquad\qquad\qquad s(0)=s_{0}\geq 0,\\\ \\\
\displaystyle\frac{de(t)}{dt}=\displaystyle{\nu s(t)i(t)}-{\rho}{e(t)},\ \ \ \
\ \ \qquad\qquad\qquad\qquad e(0)=e_{0}\geq 0,\\\ \\\
\displaystyle\frac{di(t)}{dt}=\displaystyle{\rho
e(t)}-{\delta}{i(t)}-\eta_{1}(t)i(t),\ \ \qquad\qquad\qquad i(0)=i_{0}\geq
0,\\\ \\\ \displaystyle\frac{dr(t)}{dt}=\displaystyle{\delta
i(t)+\eta_{1}(t)i(t)+\eta_{2}(t)s(t)},\ \ \ \ \qquad\qquad r(0)=r_{0}\geq
0.\end{cases}$ (35)

where the control set is given as:

${\mathcal{G}}_{5}=\\{(\eta_{1},\eta_{2})\mathrel{\mathop{\mathchar
58\relax}}\eta_{i}(t)\ \textrm{is lebesgue measurable},\ 0\leq\eta_{i}(t)\leq
1,\ t\in[0,t_{f}]\ \textrm{for}\ i=1,2\\}$

where $D_{1}$ is the weight constant on infected individuals and $D_{2}$ and
$D_{3}$ are also positive weight parameters associated with controls
$\eta_{1}(t)$ and $\eta_{2}(t)$ respectively.

The Hamiltonian function $H$ is given as

$\begin{split}H&=D_{1}i(t)]+\dfrac{D_{2}}{2}\eta^{2}_{1}(t)+\dfrac{D_{3}}{2}\eta^{2}_{2}(t)\\\
\\\ &+\varphi_{1}\left[-\nu s(t)i(t)-\eta_{2}(t)s(t)\right]\\\ \\\
&+\varphi_{2}\left[\nu s(t)i(t)-\rho e(t)\right]\\\ \\\
&+\varphi_{3}\left[\rho e(t)-\delta i(t)-\eta_{1}(t)i(t)\right]\\\ \\\
&+\varphi_{4}\left[\delta
i(t)+\eta_{1}(t)i(t)+\eta_{2}(t)s(t)\right]\end{split}$ (36)

We then apply the Maximum Principle proposed by authors in (Pontryagin et al.,
1962) to determine an optimal solution as

Assume that $(w,\eta)$ is an optimal solution for a controlled dynamical
system, then there exist adjoint vector function $\varphi=(\varphi_{1},\
\varphi_{2},\ \cdots\varphi_{n})$ which satisfy the following equations

$\begin{cases}\frac{dw}{dt}={\frac{\partial
H(t,w,\eta,\varphi)}{\partial\varphi}},\\\ \\\ 0={\frac{\partial
H(t,w,\eta,\varphi)}{\partial\eta}},\\\ \\\
\frac{d\varphi}{dt}=-{\frac{\partial H(t,w,\eta,\varphi)}{\partial
w}}.\end{cases}$ (37)

Using the formulated Hamiltonian function (36) and equation (37), the adjoint
equations and the control characterisation are presented in the following
Theorem.

###### Theorem 5

Let $\eta^{*}_{1}$ and $\eta^{*}_{2}$ be optimal control pair and optimal
state solutions $s^{*},e^{*},i^{*},r^{*}$ of the corresponding controlled
dynamical system (34)-(35) that minimize $J(\eta_{1},\eta_{2})$ over
${\mathcal{G}}_{5}$. Then there exist adjoint variables $\varphi_{i}$ for
$i=1,2,3,4$ satisfying

$\begin{cases}\displaystyle\frac{d\varphi_{1}}{dt}=\displaystyle{\nu
i^{*}(t)(\varphi_{1}-\varphi_{2})+\eta^{*}_{2}(t)(\varphi_{1}-\varphi_{4})},\qquad\\\
\\\
\displaystyle\frac{d\varphi_{2}}{dt}=\displaystyle{\rho(\varphi_{2}-\varphi_{3})},\qquad\\\
\\\ \displaystyle\frac{d\varphi_{3}}{dt}=\displaystyle{-D_{1}+\nu
s^{*}(t)(\varphi_{1}-\varphi_{2})+\delta(\varphi_{3}-\varphi_{4})+\eta^{*}_{1}(t)(\varphi_{3}-\varphi_{4})},\qquad\\\
\\\ \displaystyle\frac{d\varphi_{4}}{dt}=\displaystyle{0}.\qquad\end{cases}$
(38)

with transversality conditions $\varphi_{1}(t_{f})=0,\ \varphi_{2}(t_{f})=0,\
\varphi_{3}(t_{f})=0,\ \varphi_{4}(t_{f})=0$

and the control pair $\eta^{*}_{1}$ and $\eta^{*}_{2}$ satisfies the
optimality condition.

$\begin{cases}\displaystyle\eta^{*}_{1}(t)=min\left\\{max\left\\{0,\dfrac{i^{*}(t)}{D_{2}}(\varphi_{3}-\varphi_{4})\right\\},1\right\\}\\\
\\\
\displaystyle\eta^{*}_{2}(t)=min\left\\{max\left\\{0,\dfrac{s^{*}(t)}{D_{3}}(\varphi_{1}-\varphi_{4})\right\\},1\right\\}\end{cases}$
(39)

###  3.4 Numerical Simulations and Discussions

As we did in subsection 2.3, we have used the forward-backward sweep numerical
scheme with fourth order Runge-Kutta method to solve the optimality systems
formulated for the three optimal control strategies. For our numerical
simulations, we have adapted the same model parameter values ($\nu=0.2,\
\delta=0.1,\ \rho=0.1887$) as in the work by the authors in (Rachah and
Torres, 2016) with initial conditions: $s_{0}=0.88,\ e_{0}=0.07,\ i_{0}=0.05,\
r_{0}=0$. We have assumed $D=5,\ K_{1}=1,\ K_{2}=5,\ K_{3}=5,\ D_{1}=1,\
D_{2}=5$ and $D_{3}=5$. Figure 6 shows solution paths for susceptible
individuals with control strategies $1,2,3$ and without control. It is clear
from the plot that, there is a significant decrease in the number of
susceptible individuals with control strategies than without control. In Fig.
7, there is a rapid decrease in exposed individuals with control strategies
$1,2,3$ than without control. A similar effect can be observed in Fig. 8 for
the infected individuals. In Fig. 9, there is a rapid increase in recovered
individuals with control strategies $1,2,3$ than without control. Figures 10,
11, and 12 represents optimal control functions for strategies 1, 2 and 3
respectively.

Figure 6: Solution path for Susceptible individuals with three control
strategies and without control, $\nu=0.2$, $\rho=0.1887$, $\delta=0.1$ Figure
7: Solution paths for Exposed individuals with three control strategies and
without control, $\nu=0.2$, $\rho=0.1887$, $\delta=0.1$ Figure 8: Solution
paths for Infected individuals with three control strategies and without
control, $\nu=0.2$, $\rho=0.1887$, $\delta=0.1$ Figure 9: Solution paths for
Recovered individuals with three control strategies and without control,
$\nu=0.2$, $\rho=0.1887$, $\delta=0.1$ Figure 10: Optimal control function
$\eta(t)$ for strategy 1 Figure 11: Optimal control functions $\eta(t)$ for
strategy 2 Figure 12: Optimal control functions $\eta_{1}(t)$ and
$\eta_{2}(t)$ for strategy 3

##  4 Conclusion

In this paper, we have studied epidemiological models for Ebola virus disease
using nonlinear system of ordinary differential equation and optimal control
theory. We have used indirect methods in optimal control to study existing
mathematical models proposed by the authors in (Rachah and Torres, 2015,
2016). Using the Pontryagin’s maximum principle and Hamiltonian function, we
have derived adjoint equations and optimality system for the various optimal
control strategies. From the simulations results, we observed that, SIR model
with optimal control strategies shows significant decrease in the proportions
of infected and susceptible individuals and a rapid increase in the recovered
individuals compared to SIR model without control. A similar effect was
observed in the SEIR model with control strategies where there was a
significant decrease in Susceptible, Exposed and Infected individuals and a
rapid increase in the Recovered individuals. Our numerical results are similar
to the once generated by the authors in (Rachah and Torres, 2015, 2016) who
applied direct methods in optimal control for their mathematical models.
Therefore, following the numerical results, we can conclude that, effective
educational campaigns and vaccination of susceptible individuals as were as
effective treatments of infected individuals can help reduced the disease
transmission. In this work, we have also analytically and numerical studied
one additional optimal control strategy for the SIR epidemic model.

## References

  * Agusto et al. (2015) F. B . Agusto, M. I. Teboh-Ewungkem, and A. B Gumel.  Mathematical assessment of the effect of traditional beliefs and customs on the transmission dynamics of the 2014 ebola outbreaks.  _BMC medicine_ , 13(1):96, 2015.
  * Ahmad et al. (2016) M. D. Ahmad, M. Usman, A. Khan, and M. Imran.  Optimal control analysis of ebola disease with control strategies of quarantine and vaccination.  _Infectious diseases of poverty_ , 5(1):72, 2016\.
  * Althaus (2014) C. L. Althaus.  Estimating the reproduction number of ebola virus (ebov) during the 2014 outbreak in west africa.  _PLoS currents_ , 6, 2014.
  * Area et al. (2018) I. Area, F. Ndaïrou, J. J Nieto, C. J. Silva, and D. F. M. Torres.  Ebola model and optimal control with vaccination constraints.  _Journal of Industrial and Management Optimization_ , 14(2):427–466, 2018.
  * Atangana and Goufo (2014) A. Atangana and E. F. D. Goufo.  On the mathematical analysis of ebola hemorrhagic fever: deathly infection disease in west african countries.  _BioMed research international_ , 2014, Article ID 261383:7 pages, 2014.
  * Berge et al. (2017) T. Berge, J. M-S. Lubuma, G. M Moremedi, N. Morris, and R. Kondera-Shava.  A simple mathematical model for ebola in africa.  _Journal of biological dynamics_ , 11(1):42–74, 2017.
  * Bonyah et al. (2016) E. Bonyah, K. Badu, and S. K. Asiedu-Addo.  Optimal control application to an ebola model.  _Asian Pacific Journal of Tropical Biomedicine_ , 6(4):283–289, 2016.
  * Bonyah et al. (2017) E. Bonyah, M. A. Khan, K. O. Okosun, and S. Islam.  A theoretical model for zika virus transmission.  _PloS one_ , 12(10):e0185540, 2017.
  * Bonyah et al. (2019) E. Bonyah, M. A. Khan, K. O. Okosun, and J. F. Gómez-Aguilar.  On the co‐infection of dengue fever and zika virus.  _Optim Control Appl Meth_ , pages 1–28, 2019.
  * Chowell and Nishiura (2014) G. Chowell and H. Nishiura.  Transmission dynamics and control of ebola virus disease (evd): a review.  _BMC medicine_ , 12(1):196, 2014.
  * Chowell et al. (2004) G. Chowell, N. W. Hengartner, C. Castillo-Chavez, P. W. Fenimore, and J. M. Hyman.  The basic reproductive number of ebola and the effects of public health measures: the cases of congo and uganda.  _Journal of theoretical biology_ , 229(1):119–126, 2004.
  * Dénes and Gumel (2019) A. Dénes and A. B. Gumel.  Modeling the impact of quarantine during an outbreak of ebola virus disease.  _Infectious Disease Modelling_ , 4:12–27, 2019.
  * Diaz et al. (2018) P. Diaz, P. Constantine, K. Kalmbach, E. Jones, and S. Pankavich.  A modified seir model for the spread of ebola in western africa and metrics for resource allocation.  _Applied Mathematics and Computation_ , 324:141–155, 2018\.
  * Dong et al. (2015) F. Dong, D. Xu, Z. Wang, and M. Dong.  Evaluation of ebola spreading in west africa and decision of optimal medicine delivery strategies based on mathematical models.  _Infection, Genetics and Evolution_ , 36:35–40, 2015.
  * Goufo et al. (2016) E. F. D. Goufo, M. K. Pene, and S. Mugisha.  Stability analysis of epidemic models of ebola hemorrhagic fever with non-linear transmission.  _Journal of Nonlinear Sciences & Applications (JNSA)_, 9(6), 2016.
  * Grigorieva et al. (2017) E. V. Grigorieva, P. B. Deignan, and E. N Khailov.  Optimal control problem for a seir type model of ebola epidemics.  _Revista de Matemática: Teoría y Aplicaciones_ , 24(1):79–96, 2017.
  * Hethcote (2000) H. W. Hethcote.  The mathematics of infectious diseases.  _SIAM review_ , 42(4):599–653, 2000.
  * Jiang et al. (2017) S. Jiang, K. Wang, C. Li, G. Hong, X. Zhang, M. Shan, H. Li, and J. Wang.  Mathematical models for devising the optimal ebola virus disease eradication.  _Journal of translational medicine_ , 15(1):124, 2017.
  * Khan et al. (2015) A. Khan, M. Naveed, M. Dur-e-Ahmad, and M. Imran.  Estimating the basic reproductive ratio for the ebola outbreak in liberia and sierra leone.  _Infectious diseases of poverty_ , 4(1):13, 2015\.
  * Khan et al. (2019) M. A. Khan, S. W. Shah, S. Ullah, and J. F. Gómez-Aguilar.  A dynamical model of asymptomatic carrier zika virus with optimal control strategies.  _Nonlinear Analysis: Real World Applications_ , 50:144–170, 2019.
  * Lenhart and Workman (2007) S. Lenhart and J. T. Workman.  _Optimal control applied to biological models_.  Chapman and Hall/CRC, 2007.
  * Luo et al. (2019) D. Luo, R. Zheng, D. Wang, X. Zhang, Y. Yin, K. Wang, and W. Wang.  Effect of sexual transmission on the west africa ebola outbreak in 2014: a mathematical modelling study.  _Scientific reports_ , 9, 2019.
  * Miyaoka et al. (2019) T. Y. Miyaoka, S. Lenhart, and J. C. A. Meyer.  Optimal control of vaccination ina vector-borne reaction-diffusion model applied to zika virus.  _Journal of Mathematical Biology_ , pages 1–28, 2019.
  * Ndanguza et al. (2017) D. Ndanguza, I. S. Mbalawata, H. Haario, and J. M. Tchuenche.  Analysis of bias in an ebola epidemic model by extended kalman filter approach.  _Mathematics and Computers in Simulation_ , 142:113–129, 2017.
  * Ngwa and Teboh-Ewungkem (2016) G. A. Ngwa and M. I. Teboh-Ewungkem.  A mathematical model with quarantine states for the dynamics of ebola virus disease in human populations.  _Computational and mathematical methods in medicine_ , vol. 2016, Article ID 9352725:29 pages, 2016.
  * Pontryagin et al. (1962) L. S. Pontryagin, V.G. Boltyanskii, R. V. Gamkrelidze, and E. F. Mishchenko.  The mathematical theory of optimal processes.  _New York, John Wiley & Sons_, 1962.
  * Rachah and Torres (2015) A. Rachah and D. F. M. Torres.  Mathematical modelling, simulation, and optimal control of the 2014 ebola outbreak in west africa.  _Discrete Dynamics in Nature and Society_ , vol. 2015, Article ID 842792:9 pages, 2015.
  * Rachah and Torres (2016) A. Rachah and D. F. M. Torres.  Dynamics and optimal control of ebola transmission.  _Mathematics in Computer Science_ , 10(3):331–342, 2016.
  * Sharomi and Malik (2017) O. Sharomi and T. Malik.  Optimal control in epidemiology.  _Annals of Operations Research_ , 251(1-2):55–71, 2017.
  * Tulu et al. (2017) T. W. Tulu, B. Tian, and Z. Wu.  Mathematical modeling, analysis and markov chain monte carlo simulation of ebola epidemics.  _Results in physics_ , 7:962–968, 2017.
  * Weitz and Dushoff (2015) J. S. Weitz and J. Dushoff.  Modeling post-death transmission of ebola: challenges for inference and opportunities for control.  _Scientific reports_ , 5:8751, 2015.
  * WHO (2019) WHO.  Ebola virus disease, https://www.who.int/news-room/fact-sheets/detail/ebola-virus-disease, 2019.
  * Xia et al. (2015) Zhi-Qiang Xia, Shi-Fu Wang, Shen-Long Li, Liu-Yu Huang, Wen-Yi Zhang, Gui-Quan Sun, Zhong-Tao Gai, and Zhen Jin.  Modeling the transmission dynamics of ebola virus disease in liberia.  _Scientific reports_ , 5:13857, 2015.
  * Zakary et al. (2017) O. Zakary, M. Rachik, and I. Elmouki.  A multi-regional epidemic model for controlling the spread of ebola: awareness, treatment, and travel-blocking optimal control approaches.  _Mathematical Methods in the Applied Sciences_ , 40(4):1265–1279, 2017.