# Phonon-assisted tunneling through quantum dot systems connected to Majorana
bound states

Levente Máthé<EMAIL_ADDRESS>National Institute for R & D of
Isotopic and Molecular Technologies, 67-103 Donat, 400293 Cluj-Napoca, Romania
Faculty of Physics, Babeș-Bolyai University, 1 Kogălniceanu, 400084 Cluj-
Napoca, Romania  Zoltán Kovács-Krausz  Department of Physics, Institute of
Physics, Budapest University of Technology and Economics, Műegyetem rkp. 3.,
H-1111 Budapest, Hungary  Ioan Botiz  Faculty of Physics, Babeș-Bolyai
University, 1 Kogălniceanu, 400084 Cluj-Napoca, Romania  Interdisciplinary
Research Institute in Bio-Nano-Sciences, Babeș-Bolyai University, 42 Treboniu
Laurian, 400271 Cluj-Napoca, Romania  Ioan Grosu  Faculty of Physics, Babeș-
Bolyai University, 1 Kogălniceanu, 400084 Cluj-Napoca, Romania  Khadija El
Anouz  Laboratory of R & D in Engineering Sciences, Faculty of Sciences and
Techniques Al‑Hoceima, Abdelmalek Essaadi University, Tetouan 93000, Morocco
Abderrahim El Allati  Laboratory of R & D in Engineering Sciences, Faculty of
Sciences and Techniques Al‑Hoceima, Abdelmalek Essaadi University, Tetouan
93000, Morocco  Liviu P. Zârbo  National Institute for R & D of Isotopic and
Molecular Technologies, 67-103 Donat, 400293 Cluj-Napoca, Romania

###### Abstract

We theoretically analyze phonon-assisted tunneling transport in a quantum dot
side-connected to a Majorana bound state in a topological superconducting
nanowire. We investigate the behavior of the current through the dot, for a
range of experimentally relevant parameters, in the presence of one long-wave
optical phonon. We consider the current-gate voltage, the current-bias
voltage, and the current-dot-Majorana coupling characteristics under the
influence of the electron-phonon coupling. In the absence of electron-phonon
interaction, the Majorana bound states suppress the current when the gate
voltage matches the Fermi level, but the increase of the bias voltage
counteracts this effect. In the presence of electron-phonon coupling, the
current behaves similarly as a function of the renormalized gate voltage. As
an added feature at large bias voltages, it presents a dip or a plateau,
depending on the size of dot-Majorana coupling. Lastly, we show that the
currents are most sensitive to, and depend non-trivially on the parameters of
the Majorana circuit element, in the regime of low temperatures combined with
low voltages. Our results provide insights into the complex physics of quantum
dot devices used to probe Majorana bound states.

quantum dot; Majorana bound states; electron-phonon interaction; quantum
transport

##  I INTRODUCTION

Majorana bound states (MBSs) are zero-energy excitations in topological
materials known to form a potential platform for solid state quantum
computation due to their non-Abelian statistics [1, 2]. Previous theoretical
works [3, 4] considered devices based on semiconducting nanowires realized
from InAs or InSb with strong spin-orbit coupling located in the proximity of
$s$-wave superconductors (SCs). The latter were threaded by an external
magnetic field in order to drive the nanowire into its topological
superconducting phase by creating MBSs. Such theoretical proposals were
further confirmed by experiments [5]. Other theoretical works proposed
experimental setups realized from topological insulators [6], magnetic
nanoparticles on SCs [7], nanomagnets [8] and $p$-wave SCs [9, 10] to create
MBSs.

A minimal setup to probe MBSs in topological superconducting nanowires (TSNWs)
requires the coupling of the nanowire to a quantum dot (QD), which introduces
regular fermionic degrees of freedom [11, 12]. The presence of MBSs requires
the conductance to take the value of $e^{2}/2h$, which is measured through the
QD via normal leads [12]. Several theoretical designs based on QDs-TSNWs
setups have been considered in order to probe the MBSs via transport
properties such as (thermal) [13, 14, 15, 16, 17] conductance [18, 19, 20, 21,
22, 23, 24, 25, 26, 27, 28, 29, 30], current noise [18, 31, 32, 33, 34] and
Josephson current [35, 36]. Details on the experimental detection of MBSs via
transport characteristics measurements have been reported elsewhere [5, 37,
38, 39, 40, 41, 42]. The photon-assisted transport properties of QD-MBS setups
have been studied in the literature both theoretically [43, 44, 45] and
experimentally [46].

Over the last few years, the effect of optical phonons in transport properties
of QD-MBS systems attracted great attention [47, 48, 49, 50, 51, 52]. The
phonon-assisted transport properties of QDs coupled to MBSs have been studied
in Refs. [50, 51] in order to establish the connection between the electrical
current and heat generation in such systems. In a recent work, we analyzed the
phonon-assisted transport properties in a QD connected to a Majorana ring
structure [52]. We found the periodicity of zero-temperature linear
conductance, as a function of threading magnetic flux phase, to be independent
of the electron-phonon interaction (EPI), as well as of changes in QD energy
and finite values of the QD-MBS couplings when the Majorana wave-functions do
not overlap.

Figure 1: Schematic representation of a QD connected to an MBS located at one
of the ends of a TSNW. Here, $\eta_{1}$ and $\eta_{2}$ denote the Majorana
operators corresponding to the two MBSs in the TSNW. The dot is coupled to two
normal leads with coupling strength $\Gamma$. The electron in QD interacts
with a single long-wave optical phonon mode of frequency $\omega_{0}$. The
notations $\lambda$ and $\beta$ represent the QD-MBS and electron-phonon
coupling strengths, respectively.

In this work, we study the phonon-assisted transport properties of a QD
coupled to one MBS located at one of the ends of a TSNW. We measure the
current through the QD via metallic leads. The QD interacts with a single
long-wave optical phonon mode. We treat the EPI within a canonical
transformation which leads to the renormalization of QD energy, dot-leads and
dot-MBS couplings. Here, we consider the effect of EPI strength and
temperature on the QD-leads and QD-MBS couplings, which is usually neglected
in the literature [50, 51]. Under such conditions, we discuss in detail the
transport characteristics of the considered system for unhybridized and
hybridized MBSs.

The article is organized as follows. In Sec. II, we present the theoretical
model used by us, perform a canonical transformation on the system Hamiltonian
to eliminate the EPI. We then calculate the tunneling current using the
nonequilibrium Green’s function technique. We show and discuss the results in
Sec. III. Finally, we present our conclusions and discuss the significance of
our main results in Sec. IV.

##  II Theory

We consider a QD connected to two normal leads and to one of the ends of a
TSNW via an MBS, as is shown in Fig. 1. The normal leads allow the measuring
of current through the QD. The dot energy is tuned by the gate voltage $V_{g}$
applied to the gate electrode. The interaction between the QD electron and the
single long-wave optical phonon mode leads to phonon-assisted transport. The
Zeeman energy, $V_{Z}$, the largest energy scale in the system, is created by
the applied magnetic field which drives the nanowire into the topological
superconducting phase when the relation $V_{Z}>\sqrt{\Delta^{2}+\mu^{2}}$ is
fulfilled. Here, $\Delta$ and $\mu$ represent the SC energy gap and nanowire
chemical potential, respectively. The considered system is described by the
Hamiltonian [12, 53, 54, 55, 52]:

$H=H_{\text{leads}}+H_{\text{MBS}}+H_{\text{ph}}+H_{\text{QD}}+H_{\text{tun}}.$
(1)

Here, the Hamiltonian $H_{\text{leads}}$ models the noninteracting electrons
in both leads,

$H_{\text{leads}}=\sum_{\gamma,\mathbf{k}}\varepsilon_{\gamma\mathbf{k}}\,c_{\gamma\mathbf{k}}^{\dagger}c_{\gamma\mathbf{k}},$
(2)

where $c_{\gamma\mathbf{k}}^{\dagger}$ ($c_{\gamma\mathbf{k}}$) is the
creation (annihilation) operator for an electron with momentum $\mathbf{k}$ in
the left ($\gamma=L$) and right ($\gamma=R$) leads. Therefore,
$\varepsilon_{\gamma\mathbf{k}}=\varepsilon_{\mathbf{k}}-\mu_{\gamma}$
represent the single-particle energies and the chemical potential
$\mu_{\gamma}$. Here, the leads are at the same temperature ($T_{\gamma}=T$).
The next term in Eq. (1), $H_{\text{MBS}}$, describes the interaction between
MBSs,

$H_{\text{MBS}}=i\varepsilon_{M}\eta_{1}\eta_{2},$ (3)

where $\varepsilon_{M}\propto e^{-L/\xi}$ is the Majorana overlap energy with
the TSNW length ($L$) and superconducting coherence length ($\xi$). Here,
$\eta_{1}$ and $\eta_{2}$ are the Majorana operators related to the two MBSs
located at the opposite ends of the TSNW. The third term in Eq. (1),
$H_{\text{ph}}$, models the longitudinal optical phonon mode,

$H_{\text{ph}}=\hbar\omega_{0}a^{\dagger}a,$ (4)

with the energy $\hbar\omega_{0}$. Here, $a^{\dagger}$ and $a$ are the phonon
creation and annihilation operators. The QD Hamiltonian, $H_{\text{QD}}$, is
given by

$H_{\text{QD}}=\varepsilon_{d}d^{\dagger}d+\beta(a+a^{\dagger})d^{\dagger}d,$
(5)

where $\varepsilon_{d}$ is the QD energy and $d^{\dagger}(d)$ is the creation
(annihilation) operator for an electron in the QD. The EPI is modeled by the
second part of Eq. (5) where $\beta$ is the electron-phonon coupling strength.
The last term in Eq. (1), $H_{\text{tun}}$, represents the tunneling
Hamiltonian,

$H_{\text{tun}}=(\lambda
d-\lambda^{*}d^{\dagger})\eta_{1}+\sum_{\gamma,\mathbf{k}}\big{(}V_{\gamma\mathbf{k}}c_{\gamma\mathbf{k}}^{\dagger}d+V^{*}_{\gamma\mathbf{k}}d^{\dagger}c_{\gamma\mathbf{k}}\big{)},$
(6)

where the first term in Eq. (6) describes the coupling of strength $\lambda$
between the QD and the MBS $\eta_{1}$ located at one of the end of TSNW. The
second component of Eq. (6) refers to the coupling between the lead $\gamma$
and the QD, and is characterized by the $V_{\gamma\mathbf{k}}$ tunneling
amplitude. For further calculations, the Majorana operators $\eta_{1}$ and
$\eta_{2}$ in Eqs. (3) and (6) will be replaced with regular fermionic
operators $\eta_{1}=(f^{\dagger}+f)/\sqrt{2}$ and
$\eta_{2}=i(f^{\dagger}-f)/\sqrt{2}$. In order to explore the transport
properties of the system, we bias the QD as $\mu_{L}=-\mu_{R}=eV/2$, while the
SC is grounded, i.e., $\mu_{S}=0$. We also restrict our calculations to the
wide-band limit [56] for a symmetrically coupled QD-lead system displaying an
electron-hole symmetry, i.e.,
$\Gamma_{\gamma}^{e}=\Gamma_{\gamma}^{h}=\Gamma_{\gamma}=\Gamma$, where
$\Gamma_{\gamma}^{e(h)}=2\pi\sum_{\mathbf{k}}|V_{\gamma\mathbf{k}}|^{2}\delta(\varepsilon\mp\varepsilon_{\gamma\mathbf{k}})$
is the coupling between the dot and the lead $\gamma$ for electrons (holes).
Furthermore, we consider the relatively weak electron-phonon coupling limit
[53], by employing the nonequilibrium Green’s function technique [56, 57, 58]
in the subgap regime $|eV|<\Delta$. In this case, at finite temperature the
current takes the form [52]:

$I=\frac{ie}{2h}\Gamma\int
d\varepsilon\big{[}f_{L}^{e}(\varepsilon)-f_{R}^{e}(\varepsilon)\big{]}[\mathbf{G}_{d}^{>}(\varepsilon)-\mathbf{G}_{d}^{<}(\varepsilon)]_{11}.$
(7)

Here, $f_{\gamma}^{e}(\varepsilon)$ represents the Fermi-Dirac distribution
function for electrons in lead $\gamma$ and
$\mathbf{G}_{d}^{<(>)}(\varepsilon)$ is the lesser (greater) Green’s function
matrix of the QD in Nambu space [59, 60]. Next, we set $k_{B}=\hbar=1$. The
corresponding lesser and greater Green’s functions, appearing in Eq. (7), are
determined by using a canonical transformation $\bar{H}=e^{S}He^{-S}$ with
$S=(\beta/\omega_{0})d^{\dagger}d(a^{\dagger}-a)$ which aims to eliminate the
electron-phonon coupling term in the Hamiltonian given by Eq. (1) [53, 54, 61,
62, 52]. Thus, the transformed Hamiltonian becomes
$\bar{H}=\bar{H}_{\text{El}}+H_{\text{ph}}$ where the electron term is

$\bar{H}_{\text{El}}=H_{\text{leads}}+H_{\text{MBS}}+\tilde{H}_{\text{QD}}+\bar{H}_{\text{tun}},$
(8)

with

$\begin{split}\bar{H}_{\text{tun}}&=\frac{1}{\sqrt{2}}(\tilde{\lambda}d-\tilde{\lambda}^{*}d^{\dagger})(f+f^{\dagger})\\\
&+\sum_{\gamma,\mathbf{k}}\big{(}\tilde{V}_{\gamma\mathbf{k}}c_{\gamma\mathbf{k}}^{\dagger}d+\tilde{V}^{*}_{\gamma\mathbf{k}}d^{\dagger}c_{\gamma\mathbf{k}}\big{)}.\end{split}$
(9)

Here, $\tilde{H}_{\text{QD}}=\tilde{\varepsilon}_{d}d^{\dagger}d$ and the
renormalized dot energy $\tilde{\varepsilon}_{d}=\varepsilon_{d}-g\omega_{0}$,
with $g$ is being equal to $g=(\beta/\omega_{0})^{2}$. Furthermore, the QD-
leads and QD-MBS couplings become renormalized as
$\tilde{V}_{\gamma\mathbf{k}}=V_{\gamma\mathbf{k}}X$ and
$\tilde{\lambda}=\lambda X$ with $X=\exp(-(\beta/\omega_{0})(a^{\dagger}-a))$.
Here, we apply the approximation $X\approx\langle
X\rangle=\exp(-g(N_{\text{ph}}+\frac{1}{2}))$ [54], where
$N_{\text{ph}}=1/(e^{\omega_{0}/T}-1)$ being the Bose-Einstein distribution
function. This approximation holds if
$V_{\gamma\mathbf{k}},\,\lambda\ll\text{min}(\beta,\Delta)$ or
$\beta\ll\text{min}(V_{\gamma\mathbf{k}},\lambda,\Delta)$ [54, 63, 47]. The
lesser and greater Green’s functions are given by

$\displaystyle\mathbf{G}_{d}^{<}(\varepsilon)$
$\displaystyle=\sum_{l=-\infty}^{\infty}\mathcal{L}_{l}\mathbf{\tilde{G}}_{d}^{<}(\varepsilon+l\omega_{0}),$
(10) $\displaystyle\mathbf{G}_{d}^{>}(\varepsilon)$
$\displaystyle=\sum_{l=-\infty}^{\infty}\mathcal{L}_{l}\mathbf{\tilde{G}}_{d}^{>}(\varepsilon-l\omega_{0}),$

where
$\mathcal{L}_{l}=e^{-g(2N_{\text{ph}}+1)}e^{l\omega_{0}/(2T)}I_{l}(2g\sqrt{N_{\text{ph}}(N_{\text{ph}}+1)})$
is the Franck-Condon factor at finite temperature [54]. Here, $I_{l}(z)$ is
the lth order modified Bessel function of the first kind. Note that
$\mathcal{L}_{l}$ becomes $\mathcal{L}_{l}=e^{-g}g^{l}/l!$ for $l\geq 0$,
while $\mathcal{L}_{l}=0$ for $l<0$ at $T=0$. The dressed lesser (greater)
Green’s function $\mathbf{\tilde{G}}_{d}^{<(>)}$ is calculated by employing
the Keldysh equation
$\mathbf{\tilde{G}}_{d}^{<(>)}=\mathbf{\tilde{G}}_{d}^{r}\mathbf{\tilde{\Sigma}}^{<(>)}\mathbf{\tilde{G}}_{d}^{a}$
with the use of corresponding lesser (greater) self-energy
$\mathbf{\tilde{\Sigma}}^{<(>)}$. Thus, the current given in relation (7), at
finite temperature, reads [52]

$\begin{split}I=\frac{e}{2h}\Gamma\tilde{\Gamma}\sum_{l=-\infty}^{\infty}\mathcal{L}_{l}\int&d\varepsilon\big{[}f_{L}^{e}(\varepsilon)-f_{R}^{e}(\varepsilon)\big{]}\big{\\{}(|\tilde{G}_{d11}^{r}(\varepsilon+l\omega_{0})|^{2}+|\tilde{G}_{d12}^{r}(\varepsilon+l\omega_{0})|^{2})\big{[}f_{L}^{e}(\varepsilon+l\omega_{0})+f_{R}^{e}(\varepsilon+l\omega_{0})\big{]}\\\
&+(|\tilde{G}_{d11}^{r}(\varepsilon-l\omega_{0})|^{2}+|\tilde{G}_{d12}^{r}(\varepsilon-l\omega_{0})|^{2})\big{[}2-f_{L}^{e}(\varepsilon-l\omega_{0})-f_{R}^{e}(\varepsilon-l\omega_{0})\big{]}\big{\\}}.\\\
\end{split}$ (11)

Here, the dressed retarded Green’s functions,
$\tilde{G}_{d11}^{r}(\varepsilon)=\langle\langle
d|d^{\dagger}\rangle\rangle_{\varepsilon}^{r}$ and
$\tilde{G}_{d12}^{r}(\varepsilon)=\langle\langle
d|d\rangle\rangle_{\varepsilon}^{r}$ are calculated by employing the equation
of motion technique [64, 52]:

$\displaystyle\tilde{G}_{d11}^{r}(\varepsilon)=\frac{\varepsilon+\tilde{\varepsilon}_{d}+i\tilde{\Gamma}-|\tilde{\lambda}|^{2}K}{(\varepsilon-\tilde{\varepsilon}_{d}+i\tilde{\Gamma})(\varepsilon+\tilde{\varepsilon}_{d}+i\tilde{\Gamma})-2(\varepsilon+i\tilde{\Gamma})|\tilde{\lambda}|^{2}K},$
(12)
$\displaystyle\tilde{G}_{d12}^{r}(\varepsilon)=\frac{-|\tilde{\lambda}|^{2}K}{(\varepsilon-\tilde{\varepsilon}_{d}+i\tilde{\Gamma})(\varepsilon+\tilde{\varepsilon}_{d}+i\tilde{\Gamma})-2(\varepsilon+i\tilde{\Gamma})|\tilde{\lambda}|^{2}K},$

with $K=\varepsilon/(\varepsilon^{2}-\varepsilon_{M}^{2})$. Note that if
$\varepsilon_{M}=0$, the retarded Green’s functions given by Eq. (12) reduce
to the results of [48]. Note that the current can be simply determined at zero
temperature by replacing the Fermi-Dirac function $f_{\gamma}(x)$ in Eq. (11)
with the Heaviside one $\theta(\mu_{\gamma}-x)$.

##  III Results and discussion

In the following, we discuss the transport properties of the QD-MBS system
introduced above for a few experimentally relevant parameter regimes.

As already mentioned in Sec. II, the system parameters must be smaller than
the SC energy gap $\Delta$, which is on the order of $250\,\mu eV$ in TSNW
experiments [5]. In addition, our phonon-assisted transport calculations are
performed in the limit where the QD-lead and QD-MBS couplings are weaker than
the the electron-phonon coupling strength [52]. In our calculations, the
symmetrical QD-lead coupling $\Gamma$ is used as the energy unit. The optical
phonon energy $\omega_{0}$ and EPI strength are considered to be
$\omega_{0}=5\,\Gamma$ and $\beta=2.5\,\Gamma$. For more details regarding of
the choice of parameters based on experimental measurements, see Ref. [52] and
references therein.

We analyze the transport characteristics of our system at a finite
temperature, thus serving as a relevant case for real systems. We also
consider the $\beta$ and $T$ dependence of $\tilde{\Gamma}=\Gamma
e^{-g(2N_{ph}+1)}$ and $|\tilde{\lambda}|=|\lambda|e^{-g(N_{ph}+1/2)}$. Recall
that the QD energy can be tuned by the gate voltage ($V_{g}$) applied to the
gate electrode, i.e., $\varepsilon_{d}\propto V_{g}$. In the following
subsections, we show how the transport currents are influenced by the system
parameters.

i) The effect of QD-MBS coupling $|\lambda|$ and bias voltage $V$ on current
vs gate voltage characteristics

We first study the effect of QD-MBS coupling $|\lambda|$ and bias voltage $eV$
on the characteristics of current vs
$\tilde{\varepsilon}_{d}=\varepsilon_{d}-\beta^{2}/\omega_{0}$ for
unhybridized Majoranas in the presence of EPI at a finite temperature
$T=0.1\,\Gamma$. Fig. 2(a) shows the current $I$ as a function of
$\tilde{\varepsilon}_{d}$ for different values of the dot-MBS coupling
$|\lambda|$ when the system is biased as $eV=2\,\Gamma$. The calculations are
made for unhybridized MBSs ($\varepsilon_{M}=0$) in absence and presence of
EPI with a fixed electron-phonon coupling strength $\beta=2.5\,\Gamma$. We
observe that in absence of EPI and MBSs, a single Lorentzian resonant peak
emerges at $\tilde{\varepsilon}_{d}=\varepsilon_{d}=0$. The amplitude of this
peak is reduced when the dot couples to one of the ends of a TSNW (i.e.,
$|\lambda|\neq 0$) in the dot energy range
$-(\Gamma+\frac{|eV|}{2})\lesssim\varepsilon_{d}\lesssim(\Gamma+\frac{|eV|}{2})$.
Beside the dot energy domain
$|\varepsilon_{d}|\gtrsim(\Gamma+\frac{|eV|}{2})$, the magnitude of the
current $|I|$ increases slightly with $|\lambda|$ [see Fig. 2(a), dotted
lines]. In the presence of EPI, the current has a maximum at
$\tilde{\varepsilon}_{d}=0$ and the absolute value of its amplitude is reduced
compared to the $\beta=0$ case. When the QD hybridizes with the MBS
($|\lambda|\neq 0$), the spectrum of $I$ is changed [see Fig. 2(a), solid
lines]. Similarly to the no EPI case, the suppression of $|I|$ is realized
when the renormalized dot energy $\tilde{\varepsilon}_{d}$ is situated within
the energy region
$-(\tilde{\Gamma}+\frac{|eV|}{2})\lesssim\tilde{\varepsilon}_{d}\lesssim(\tilde{\Gamma}+\frac{|eV|}{2})$.
Otherwise, when
$|\tilde{\varepsilon}_{d}|\gtrsim(\tilde{\Gamma}+\frac{|eV|}{2})$, a slight
increase of the current magnitude $|I|$ with $|\lambda|$ is observed. This
behavior of the current agrees qualitatively with the results of Ref. [50]
where the effect of EPI and temperature on $\Gamma$ and $|\lambda|$ is
neglected by considering the couplings as constants. The effect of the bias
voltage $eV$ on the characteristics of $I-\tilde{\varepsilon}_{d}$ is shown
without EPI in Fig. 2(b), and with EPI of strength $\beta=2.5\,\Gamma$ in Fig.
2(c). In the absence of EPI and MBSs [Fig. 2(b), black lines], the magnitude
of $|I|$ increases with the increase of bias $|eV|$, in agreement with the
literature [65]. But, when the MBS is introduced in the system with $\beta=0$
[see Fig. 2(b), red and green lines], the magnitude of the current is reduced
with respect to the case of $|\lambda|=0$ within the dot energy region
$-(\Gamma+\frac{|eV|}{2})\lesssim\varepsilon_{d}\lesssim(\Gamma+\frac{|eV|}{2})$,
in agreement with the result shown in Fig. 2(a). In the presence of EPI with
MBSs [see Fig. 2(c)], the current responds in the same way to the change in
$|\lambda|$ for small values of the voltage $|eV|$ as in the $\beta=0$ case
under the mappings $\varepsilon_{d}\to\tilde{\varepsilon}_{d}$ and
$\Gamma\to\tilde{\Gamma}$, respectively. In the $\beta\neq 0$ case, the effect
of the hybridization of MBS with QD on the $I-\tilde{\varepsilon}_{d}$ curves
alters depending on the bias voltage [discussed also in Fig. 5(a) below].
Consequently, the change in current magnitude due to the QD-Majorana coupling
$|\lambda|$ can be counteracted by increasing the bias voltage.

Figure 2: (a) The current $I$ as a function of $\tilde{\varepsilon}_{d}$ for
different values of the QD-MBS coupling $|\lambda|$ with unhybridized MBSs
($\varepsilon_{M}=0\,\Gamma$). The bias voltage is fixed as $eV=2\,\Gamma$,
while the solid (dotted) lines correspond to the $\beta=2.5\,\Gamma$
($\beta=0\,\Gamma$) case. The current $I$ as a function of
$\tilde{\varepsilon}_{d}$ for unhybridized MBSs ($\varepsilon_{M}=0\,\Gamma$)
at different values of bias voltage $eV$ and QD-MBS coupling $|\lambda|$ in
the (b) absence and (c) presence of EPI with $\beta=2.5\,\Gamma$. Here, the
solid, dashed, dot-dashed and dotted lines correspond to voltages $eV$ equal
to $0.5\,\Gamma$, $2\,\Gamma$, $4\,\Gamma$ and $8\,\Gamma$, respectively. In
all cases the temperature is fixed at $T=0.1\,\Gamma$. Figure 3: The current
$I$ as a function of $\tilde{\varepsilon}_{d}$ for different values of the
overlap energy $\varepsilon_{M}$ with fixed electron-phonon coupling strength
$\beta=2.5\,\Gamma$ at two values of the QD-MBS coupling $|\lambda|$: (a)
$|\lambda|=0.5\,\Gamma$ and (b) $|\lambda|=1\,\Gamma$. The solid (dotted)
lines correspond to the EPI presence (absence) case. The temperature and bias
voltage are $T=0.1\,\Gamma$ and $eV=2\,\Gamma$, respectively. Insets show
current $I$ as a function of overlap energy $\varepsilon_{M}$ at
$\tilde{\varepsilon}_{d}=0\,\Gamma$ in the presence of EPI with
$\beta=2.5\,\Gamma$ (solid line) and at $\varepsilon_{d}=0$ in the absence of
it (dotted line) for QD-MBS couplings: (a) $|\lambda|=0.5\,\Gamma|$ and (b)
$|\lambda|=1\,\Gamma$. Figure 4: The current $I$ as a function of QD-MBS
coupling $|\lambda|$ at different values of the overlap energy
$\varepsilon_{M}$ for bias voltage (a) $eV=2\,\Gamma$ and (b) $eV=4\,\Gamma$.
The solid (dotted) lines correspond to the $\beta=2.5\,\Gamma$ at
$\tilde{\varepsilon}_{d}=0\,\Gamma$ ($\beta=0\,\Gamma$ at
$\varepsilon_{d}=0\,\Gamma$) case. The insets in each panel zoom in on
current. (c) The current $I$ as a function of bias voltage $eV$ and QD-MBS
coupling $|\lambda|$ for unhybridized MBSs ($\varepsilon_{M}=0\,\Gamma$) in
the absence of EPI. (d) The current $I$ as a function of overlap energy
$\varepsilon_{M}$ and QD-MBS coupling $|\lambda|$ at bias voltage
$eV=2\,\Gamma$ in the absence of EPI. The QD energy is
$\varepsilon_{d}=0\,\Gamma$ in (c) and (d). The considered temperature is
$T=0.1\,\Gamma$.

ii) The effect of MBS overlap energy $\varepsilon_{M}$ on current vs gate
voltage characteristics

Next, we investigate the effect of the overlap energy $\varepsilon_{M}$ on the
$I-\tilde{\varepsilon}_{d}$ characteristics. Fig. 3 shows the results for the
current $I$ as a function of QD energy $\varepsilon_{d}$ in the absence of EPI
and as a function of $\tilde{\varepsilon}_{d}$ in the presence of EPI with
coupling $\beta=2.5\,\Gamma$ at different values of the overlap energy
$\varepsilon_{M}$ for two QD-MBS coupling $|\lambda|$ values. The temperature
and bias voltage are fixed at $T=0.1\,\Gamma$ and $eV=2\,\Gamma$,
respectively. We observe that the current magnitude around at
$\tilde{\varepsilon}_{d}=\varepsilon_{d}=0$ without EPI and at
$\tilde{\varepsilon}_{d}=0$ with EPI reduces with the increase of the overlap
energy $\varepsilon_{M}$ when the QD weakly couples to the MBS [see Fig.
3(a)]. In the weak $|\lambda|$ case with a given $\varepsilon_{M}\neq 0$, we
also see that by moving away from $\varepsilon_{d}=0$ for $\beta=0$ (or
$\tilde{\varepsilon}_{d}=0$ for $\beta\neq 0$) to larger $|\varepsilon_{d}|$
(or $|\tilde{\varepsilon}_{d}|$) values, there is a critical value above which
the current changes from a reduction to an enhancement relative to
$\varepsilon_{M}=0$ curve. Therefore, we observe that the further increase of
$\varepsilon_{M}$ leads to the current amplitude at
$\tilde{\varepsilon}_{d}=0$ with EPI (or at $\varepsilon_{d}=0$ without EPI)
begins to increase [see the inset in Fig. 3(a)]. At stronger $|\lambda|$ [see
Fig. 3(b) with its inset], the current reacts differently to the change in
$\varepsilon_{M}$ relative to the $|\lambda|=0.5\,\Gamma$ case. Namely, near
$\varepsilon_{d}=0$ for $\beta=0$ or $\tilde{\varepsilon}_{d}=0$ for
$\beta\neq 0$, the current amplitude increases with $\varepsilon_{M}$. When
the QD level passes a critical energy value, the current magnitude at a given
$\varepsilon_{M}$ will be reduced relative to the $\varepsilon_{M}=0$ case
(not shown here). This critical dot energy value moves to smaller
$|\varepsilon_{d}|$ with the increase of $\varepsilon_{M}$. Note that the
current curve peak in the presence of EPI (Fig. 3, solid lines) is narrower
than the one without EPI (Fig. 3, dotted lines) because of the renormalized
QD-lead coupling $\tilde{\Gamma}$.

To further understand the regime $eV\neq 0$, $\varepsilon_{d}\neq 0$ and
$|\lambda|\neq 0$ with $\varepsilon_{M}\neq 0$, we plot in Figs. 4(a) and (b)
the current $I$ for different values of the overlap energy $\varepsilon_{M}$
and bias voltage $eV$ in the absence of EPI at $\varepsilon_{d}=0$ and in the
presence of EPI at $\tilde{\varepsilon}_{d}=0$ with electron-phonon coupling
strength $\beta=2.5\,\Gamma$, respectively. These results are represented as a
function of QD-MBS coupling $|\lambda|$ at a finite temperature
$T=0.1\,\Gamma$. We can see that in the absence of EPI with unhybridized
Majoranas [see Figs. 4(a)-(b), blue dotted lines], the maximum of the current
which emerges at $\varepsilon_{d}=0$ significantly reduces when the dot
hybridizes with the MBS, in agreement with the result from Fig. 2(a). We can
observe that $|I|$ increases with $|eV|$ at fixed values of $|\lambda|$, in
agreement with the results from Fig. 2(b). Therefore, the value of QD-MBS
coupling $|\lambda|$ where $|I|$ presents a significant reduction, shifts to
higher values of $|\lambda|$ with the increase of bias $|eV|$. In the case of
hybridized MBSs, at a given $\varepsilon_{M}$, there is a critical value for
QD-MBS coupling, where the effect of $\varepsilon_{M}$ on the current
amplitude changes from a reduction to an enhancement relative to the
$\varepsilon_{M}=0$ case. This critical value of $|\lambda|$ depends on the
bias voltage $eV$. The increasing bias $|eV|$ leads to negligible influence of
the finite overlap energies considered here on the current curves. Note here
that a relation for the current $I(eV)=-I(-eV)$ can be established, as seen
also in Fig. 5 below. To better understand the current-bias voltage dependence
when the dot-Majorana coupling changes, we plot in Fig. 4(c) the current $I$
as a function of $eV$ and $|\lambda|$ in the absence of EPI for unhybridized
MBSs at $\varepsilon_{d}=0$ and a finite temperature $T=0.1\,\Gamma$. In Fig.
4(c), we notice that the line $|\lambda|\propto eV$ represents an inflection
point which corresponds to a peak in the differential conductance. Fig. 4(d)
shows the results for the current $I$ as a function of overlap energy
$\varepsilon_{M}$ and QD-MBS coupling strength $|\lambda|$ in the absence of
EPI at a fixed voltage $eV=2\,\Gamma$, dot energy $\varepsilon_{d}=0$ and
temperature $T=0.1\,\Gamma$. The current map details further changes with the
enhancement of the MBS-MBS coupling strength $\varepsilon_{M}$ predominantly
when the QD strongly couples to the MBS [see also Fig. 4(a)]. In the presence
of EPI [see Figs. 4(a) and (b), solid lines], the amplitude of current $|I|$
is reduced relative to the $\beta=0$ case which is more visible at voltages
near $eV\approx 4\,\Gamma$ [see also Fig. 5 below for a larger $eV$ domain].
The $I-|\lambda|$ curves for $\beta\neq 0$ show the same behavior as those for
$\beta=0$ at low bias voltages.

iii) The effect of QD-MBS coupling $|\lambda|$ and MBS overlap energy
$\varepsilon_{M}$ on current vs bias voltage characteristics

In the following, we investigate the current-bias voltage characteristics of
the proposed system at different values of the QD-MBS coupling $|\lambda|$ in
the absence and presence of EPI for both unhybridized and hybridized MBSs. The
results are shown in Fig. 5 at $\tilde{\varepsilon}_{d}=0$ with EPI of
strength $\beta=2.5\,\Gamma$ and at $\varepsilon_{d}=0$ without EPI, and at a
finite temperature $T=0.1\,\Gamma$. The current $I$ shows a step-like
structure as a function of bias voltage $eV$. This structure is explained
below. In the absence of EPI and MBSs [see Fig. 5(a), red dotted line], when
the system is positively biased ($eV>0$) and $\mu_{R}<\varepsilon_{d}<\mu_{L}$
(with $\mu_{L}=-\mu_{R}=eV/2$), the dot is able to receive an electron from
lead $L$ and transfer it to lead $R$, which results in a current passing
through the dot. Such a system (without MBS or EPI) has already been detailed
in Ref. [66]. When the dot couples to the MBS, the magnitude of the current
$|I|$ reduces in the vicinity of zero-bias voltage with the increase of QD-MBS
coupling $|\lambda|$ for unhybridized MBSs [see Fig. 5(a) with its inset,
dotted lines]. The width of this voltage window, where the current is affected
by the QD-MBS coupling, becomes larger with enhancing of $|\lambda|$, in
agreement with the findings from Fig. 4(c). In the presence of EPI and absence
of MBSs [see Fig. 5(a), black solid line], the $I-V$ curve is visibly
different and new steps show up in the spectrum which correspond to the
opening of phonon-assisted tunneling channels [67]. When the dot hybridizes
with the MBS in $\beta\neq 0$ case for $\varepsilon_{M}=0$ [see Fig. 5(a),
blue and green solid lines], further changes in the $I-V$ characteristics can
be observed. For instance, in the positive bias domain, i.e., $eV\gtrsim 0$,
we see that the amplitude of $I$ decreases with the enhancement of QD-MBS
coupling $|\lambda|$ when the voltage is approximately within
$2l\omega_{0}\lesssim eV\lesssim(2l+1)\omega_{0}$. Instead, this amplitude
increases with $|\lambda|$ in the bias regimes $(2l+1)\omega_{0}\lesssim
eV\lesssim 2(l+1)\omega_{0}$ with $l=0,1,2,..$. These findings are consistent
with the results reported in Ref. [50]. For negative voltages ($eV\lesssim
0$), the changes in current are similar, taking into account the antisymmetric
nature of the $I-V$ curves. Thus, the current $|I|$ decreases with the
increase of $|\lambda|$ when $eV$ is within $(2l+1)\omega_{0}\lesssim
eV\lesssim 2(l+1)\omega_{0}$ and it increases with $|\lambda|$ when
$2l\omega_{0}\lesssim eV\lesssim(2l+1)\omega_{0}$ with $l=-1,-2,...$,
respectively.

Figure 5: (a) The current $I$ as a function of bias voltage $eV$ at different
values of QD-MBS coupling $|\lambda|$ for unhybridized MBSs
($\varepsilon_{M}=0\,\Gamma$). The current $I$ as a function of bias voltage
$eV$ for hybridized MBSs when the overlap energy $\varepsilon_{M}$ is varied
for QD-MBS couplings: (b) $|\lambda|=0.5\,\Gamma$ and (c)
$|\lambda|=1\,\Gamma$, respectively. The solid and dotted lines correspond to
the presence and absence of EPI. The electron-phonon coupling strength and
renormalized dot energy are $\beta=2.5\,\Gamma$ and
$\tilde{\varepsilon}_{d}=0\,\Gamma$, respectively. In the absence of EPI the
QD energy is $\varepsilon_{d}=0\,\Gamma$. The temperature is $T=0.1\,\Gamma$.
The insets in each panel zoom in on current near the zero-bias at positive
voltage.

The effect of overlap energy $\varepsilon_{M}$ on current-bias voltage
characteristics is illustrated in Figs. 5(b) and (c) for two values of the
dot-MBS coupling $|\lambda|$, both in the absence at $\varepsilon_{d}=0$ and
presence of EPI with strength $\beta=2.5\,\Gamma$ at
$\tilde{\varepsilon}_{d}=0$. In the absence of EPI [see Figs. 5(b) and (c),
dotted lines], we observe that the current magnitude $|I|$ increases with the
overlap energy $\varepsilon_{M}$ near the zero-bias regime, when the voltage
is constrained within $|eV|\lesssim\mathcal{V}$. Here,
$\mathcal{V}\propto(2|\lambda|+\varepsilon_{M})$ for a given $|\lambda|\neq
0$. Otherwise ($|eV|\gtrsim\mathcal{V}$), the current $|I|$ reduces slightly
with increasing $\varepsilon_{M}$. In the presence of EPI for hybridized
Majoranas [see Figs. 5(b) and (c), solid lines], the finite $\varepsilon_{M}$
significantly influences the current-bias voltage characteristics. Namely, for
a positively biased QD system ($eV\gtrsim 0$), the magnitude of the current
$|I|$ increases with $\varepsilon_{M}$ when the voltage is located within
$2l\omega_{0}\lesssim eV\lesssim 2l\omega_{0}+\mathcal{\tilde{V}}$ and
$(2l+1)\omega_{0}\lesssim eV\lesssim 2(l+1)\omega_{0}-\mathcal{\tilde{V}}$,
and decreases with the increase of $\varepsilon_{M}$ when $eV$ is within
$2l\omega_{0}+\mathcal{\tilde{V}}\lesssim eV\lesssim(2l+1)\omega_{0}$ and
$2(l+1)\omega_{0}-\mathcal{\tilde{V}}\lesssim eV\lesssim 2(l+1)\omega_{0}$,
with $l=0,1,2,...$. Here, we introduced the notation
$\mathcal{\tilde{V}}\propto(2|\tilde{\lambda}|+\varepsilon_{M})$ for a given
$|\tilde{\lambda}|\neq 0$ with the restriction
$0<\mathcal{\tilde{V}}<\omega_{0}$. For a negatively biased QD ($eV\lesssim
0$), the alteration of the current behavior as a response to the change in
$\varepsilon_{M}$ is expressed similarly to the case of $eV\gtrsim 0$ by
taking into account the asymmetric property of $I-V$ characteristics.
Consequently, when the dot couples to the MBS, the current can be amplified or
reduced by changing the value of Majorana overlap energy at a fixed bias
voltage.

Figure 6: (a) The current $I$ as a function of temperature $T$ at different
values of the QD-MBS coupling $|\lambda|$ for unhybridized MBSs
($\varepsilon_{M}=0\,\Gamma$) at bias voltage $eV=2\,\Gamma$. The inset in (a)
shows the results at $eV=0.5\,\Gamma$. The current $I$ as a function of
temperature $T$ at different values of the overlap energy $\varepsilon_{M}$
and at bias voltage $eV=2\,\Gamma$ for QD-MBS couplings: (b)
$|\lambda|=0.5\,\Gamma$ and (c) $|\lambda|=1\,\Gamma$. The insets in (b) and
(c) show the results at bias $eV=0.5\,\Gamma$. The solid (dotted) lines
represent the results with (without) EPI. The renormalized QD energy is
$\tilde{\varepsilon}_{d}=0\,\Gamma$ in the presence of EPI while the dot
energy is $\varepsilon_{d}=0\,\Gamma$ in its absence. The electron-phonon
coupling strength is $\beta=2.5\,\Gamma$.

iv) The effect of dot-MBS coupling $|\lambda|$ and Majorana overlap energy
$\varepsilon_{M}$ on current vs temperature characteristics

We plot in Fig. 6 the current $I$ as a function of temperature $T$ for
different values of the bias $eV$, QD-MBS coupling $|\lambda|$ and overlap
energy $\varepsilon_{M}$ in the absence at $\varepsilon_{d}=0$ and in the
presence of EPI with $\beta=2.5\,\Gamma$ at $\tilde{\varepsilon}_{d}=0$. We
first observe that the current is suppressed for $\beta\neq 0$ relative to its
value at $\beta=0$ [see Figs. 6(a)-(c) and their inset, solid and dotted
lines]. We find that the current shows a nonmonotonic behavior with
temperature $T$ for unhybridized Majoranas when the QD weakly couples to the
MBS at low bias voltage in both $\beta=0$ and $\beta\neq 0$ cases [inset in
Fig. 6(a), green solid and dotted lines]. Namely, $|I|$ increases with $T$ up
to a specific value, and begins decreasing above that value. With the increase
of $|\lambda|$ this tendency of the $I-T$ characteristics starts to vanish
[inset in Fig. 6(a), brown solid and dotted lines]. As the bias increases, the
nonmonotonic behavior of the current with $T$ emerges at stronger coupling
$|\lambda|$ [see Fig. 6(a), brown solid and dotted lines]. Note here that this
tendency of the $I-T$ characteristics is disappeared at higher voltages for
values of $|\lambda|$ which satisfy the approximation applied here. In the
case of hybridized MBSs with relatively strong overlap energy, the
nonmonotonic behavior of the current with temperature disappears at small
biases $eV$ [see the inset in Fig. 6(b), red and green lines]. In addition,
$|I|$ varies nonmonotonically with temperature at weaker $|\lambda|$ with
stronger $\varepsilon_{M}$ at bias $eV=2\,\Gamma$ [Fig. 6(b), green solid and
dotted lines]. This behavior of the $I-T$ curves starts to vanish with further
increase of $\varepsilon_{M}$ (not shown here). When the QD couples strongly
to the MBS, the current decreases with temperature at low bias voltages
[$eV=0.5\,\Gamma$, see the inset in Fig. 6(c)]. The nonmonotonic behavior of
the current-temperature curves holds for values of $\varepsilon_{M}$
considered here at $eV=2\,\Gamma$ and $|\lambda|=1\,\Gamma$ [see Fig. 6(c)].
Similarly to the $|\lambda|=0.5\,\Gamma$ case, the further increase of
$\varepsilon_{M}$ smears the nonmonotonic behavior of the current as a
function of temperature at $eV=2\,\Gamma$ (not shown here). Note here that the
nonmonotonic behavior of current with temperature vanishes at higher bias
voltages regardless of the $|\lambda|$ coupling values within the limits set
by the approximations used here. In this case, the overlap energy also shows
less influence on the $I-T$ characteristics. Consequently, the response of the
current to the changes in temperature is altered nontrivially depending on
bias voltage, QD-MBS coupling and Majorana overlap energy.

##  IV Summary

In the present work, we investigated the complex physics of a QD coupled to a
Majorana bound state located at one of the edges of a TSNW. In addition, the
phonon-assisted transport properties of the considered setup were explored
when the localized electrons in the QD interact with a single long-wave
optical phonon mode. When determining the current, the electron-phonon
interaction was treated by employing a canonical transformation within the
nonequilibrium Green’s function formalism. We discussed in detail the effect
of EPI on current vs gate voltage, current vs bias voltage and current vs dot-
Majorana coupling characteristics for unhybridized and hybridized MBSs at
finite temperature. We established that in the absence of EPI the dot-Majorana
coupling strength suppresses the current when the dot energy locates near the
Fermi level, especially at low bias voltages, therefore the Majorana overlap
energy and dot-Majorana coupling have a more significant impact on transport
under this low bias regime. The effect can be counteracted by increasing the
bias voltage. In the presence of EPI, the effect of dot-MBS coupling on the
current-gate voltage characteristics can be regulated by changing the bias
voltage. The current-bias voltage curves present a step-like structure in the
presence of electron-phonon coupling due to the phonon-assisted tunneling
through the dot. The effect of Majorana overlap energy on current vs bias
voltage characteristics alters depending on the bias voltage value. Note that
the current is insensitive to charge fluctuations at high voltages. We found
that the current shows a nonmonotonic behavior with temperature depending on
the values of QD-MBS coupling, overlap energy, gate and bias voltages. The
current shows sizable changes at low temperatures if the bias voltage is low.
Note that the current-temperature dependence is strongly affected even by
lower gate voltage variations. In the future, we plan to extend this
investigation to determine if such systems possess parameter regimes for which
it is easy to establish the presence or the absence of Majorana bound states
in the system via transport measurements. Finally, the device geometry
investigated in this work should be experimentally realizable by taking into
consideration the recent advancements in the field [41, 68]. Our investigation
can serve as a guide for experiments probing MBSs with QDs, helping enlarge
the understanding of topological quantum computation.

###### Acknowledgements.

The authors would like to thank Dr. Doru Sticleţ, Dr. Luiza Buimaga-Iarinca,
Larisa-Milena Pioraș-Țimbolmaș and Pál-Attila Máthé for valuable discussions.
L. M. and L. P. Z. acknowledge financial support from the MCID through the
”Nucleu” Program within the National Plan for Research, Development and
Innovation 2022-2027, project PN 23 24 01 04 and through Program 1 -
Development of the National Research and Development System, Subprogram 1.2 -
Institutional Performance - Funding Projects for Excellence in RDI, Contract
No. 37PFE/30.12.2021.

## References

  * Kitaev [2001] A. Y. Kitaev, Phys. Usp. 44, 131 (2001).
  * Kitaev [2006] A. Kitaev, Ann. Phys. 321, 2 (2006).
  * Lutchyn  _et al._ [2010] R. M. Lutchyn, J. D. Sau, and S. Das Sarma, Phys. Rev. Lett. 105, 077001 (2010).
  * Oreg  _et al._ [2010] Y. Oreg, G. Refael, and F. von Oppen, Phys. Rev. Lett. 105, 177002 (2010).
  * Mourik  _et al._ [2012] V. Mourik, K. Zuo, S. M. Frolov, S. R. Plissard, E. P. A. M. Bakkers, and L. P. Kouwenhoven, Science 336, 1003 (2012).
  * Fu and Kane [2008] L. Fu and C. L. Kane, Phys. Rev. Lett. 100, 096407 (2008).
  * Choy  _et al._ [2011] T.-P. Choy, J. M. Edge, A. R. Akhmerov, and C. W. J. Beenakker, Phys. Rev. B 84, 195442 (2011).
  * Kjaergaard  _et al._ [2012] M. Kjaergaard, K. Wölms, and K. Flensberg, Phys. Rev. B 85, 020503 (2012).
  * Read and Green [2000] N. Read and D. Green, Phys. Rev. B 61, 10267 (2000).
  * Stone and Roy [2004] M. Stone and R. Roy, Phys. Rev. B 69, 184511 (2004).
  * Flensberg [2011] K. Flensberg, Phys. Rev. Lett. 106, 090503 (2011).
  * Liu and Baranger [2011] D. E. Liu and H. U. Baranger, Phys. Rev. B 84, 201308 (2011).
  * Leijnse [2014] M. Leijnse, New J. Phys. 16, 015029 (2014).
  * López  _et al._ [2014] R. López, M. Lee, L. Serra, and J. S. Lim, Phys. Rev. B 89, 205418 (2014).
  * Khim  _et al._ [2015] H. Khim, R. López, J. S. Lim, and M. Lee, Eur. Phys. J. B 88, 015029 (2015).
  * Chi  _et al._ [2020a] F. Chi, Z.-G. Fu, J. Liu, K.-M. Li, Z. Wang, and P. Zhang, Nanoscale Res. Lett. 15, 79 (2020a).
  * Grosu and Tifrea [2023] I. Grosu and I. Tifrea, Physica E 145, 115503 (2023).
  * Cao  _et al._ [2012] Y. Cao, P. Wang, G. Xiong, M. Gong, and X.-Q. Li, Phys. Rev. B 86, 115311 (2012).
  * Lee  _et al._ [2013] M. Lee, J. S. Lim, and R. López, Phys. Rev. B 87, 241402 (2013).
  * Vernek  _et al._ [2014] E. Vernek, P. H. Penteado, A. C. Seridonio, and J. C. Egues, Phys. Rev. B 89, 165314 (2014).
  * Dessotti  _et al._ [2014] F. A. Dessotti, L. S. Ricco, M. de Souza, F. M. Souza, and A. C. Seridonio, J. Appl. Phys. 116, 173701 (2014).
  * Stefański [2015] P. Stefański, Acta Phys. Pol. A 127, 198 (2015).
  * Zeng  _et al._ [2016] Q.-B. Zeng, S. Chen, L. You, and R. Lü, Front. Phys. 12, 127302 (2016).
  * Calle  _et al._ [2020] A. M. Calle, M. Pacheco, P. A. Orellana, and J. A. Otálora, Ann. Phys. 532, 1900409 (2020).
  * Ricco  _et al._ [2020] L. S. Ricco, Y. Marques, J. E. Sanches, I. A. Shelykh, and A. C. Seridonio, Phys. Rev. B 102, 165104 (2020).
  * Majek and Weymann [2021] P. Majek and I. Weymann, Phys. Rev. B 104, 085416 (2021).
  * Smirnov [2021] S. Smirnov, Phys. Rev. B 104, 205406 (2021).
  * Gong  _et al._ [2022a] T. Gong, L.-L. Zhang, X.-F. Dai, C. Jiang, and W.-J. Gong, Eur. Phys. J. Plus 137, 122 (2022a).
  * Gong  _et al._ [2022b] T. Gong, L.-L. Zhang, S.-F. Zhang, C. Jiang, and W.-J. Gong, Ann. Phys. 447, 169143 (2022b).
  * Jiang  _et al._ [2022] C. Jiang, G.-Y. Meng, T. Gong, and L.-L. Zhang, Appl. Phys. A 128, 533 (2022).
  * Lü  _et al._ [2012] H.-F. Lü, H.-Z. Lu, and S.-Q. Shen, Phys. Rev. B 86, 075318 (2012).
  * Chen  _et al._ [2014] Q. Chen, K.-Q. Chen, and H.-K. Zhao, J. Phys.: Condens. Matter 26, 315011 (2014).
  * Lü  _et al._ [2014] H.-F. Lü, H.-Z. Lu, and S.-Q. Shen, Phys. Rev. B 90, 195404 (2014).
  * Smirnov [2022] S. Smirnov, Phys. Rev. B 105, 205430 (2022).
  * Lee and Lee [2016] Y.-L. Lee and Y.-W. Lee, Phys. Rev. B 93, 184502 (2016).
  * Chi  _et al._ [2023] F. Chi, Q.-S. Jia, J. Liu, Q.-G. Gao, Z.-C. Yi, and L.-M. Liu, Coatings 13, 10.3390/coatings13030612 (2023).
  * Finck  _et al._ [2013] A. D. K. Finck, D. J. Van Harlingen, P. K. Mohseni, K. Jung, and X. Li, Phys. Rev. Lett. 110, 126406 (2013).
  * Chen  _et al._ [2017] J. Chen, P. Yu, J. Stenger, M. Hocevar, D. Car, S. R. Plissard, E. P. A. M. Bakkers, T. D. Stanescu, and S. M. Frolov, Sci. Adv. 3, e1701476 (2017).
  * Zhang  _et al._ [2017] H. Zhang, O. Gül, S. Conesa-Boj, M. P. Nowak, M. Wimmer, K. Zuo, V. Mourik, F. K. de Vries, J. van Veen, M. W. A. de Moor, J. D. S. Bommer, D. J. van Woerkom, D. Car, S. R. Plissard, E. P. A. M. Bakkers, M. Quintero-Pérez, M. C. Cassidy, S. Koelling, S. Goswami, K. Watanabe, T. Taniguchi, and L. P. Kouwenhoven, Nat. Commun. 8, 16025 (2017).
  * Lutchyn  _et al._ [2018] R. M. Lutchyn, E. P. A. M. Bakkers, L. P. Kouwenhoven, P. Krogstrup, C. M. Marcus, and Y. Oreg, Nat. Rev. Mater. 3, 52 (2018).
  * Deng  _et al._ [2016] M. T. Deng, S. Vaitiekenas, E. B. Hansen, J. Danon, M. Leijnse, K. Flensberg, J. Nygård, P. Krogstrup, and C. M. Marcus, Science 354, 1557 (2016).
  * Sherman  _et al._ [2017] D. Sherman, J. S. Yodh, S. M. Albrecht, J. Nygård, P. Krogstrup, and C. M. Marcus, Nat. Nanotechnol. 12, 212 (2017).
  * Tang  _et al._ [2015] H.-Z. Tang, Y.-T. Zhang, and J.-J. Liu, AIP Advances 5, 127129 (2015).
  * Chi  _et al._ [2020b] F. Chi, T.-Y. He, J. Wang, Z.-G. Fu, L.-M. Liu, P. Liu, and P. Zhang, Front. Phys. (Lausanne) 8, 254 (2020b).
  * Ricco  _et al._ [2022] L. S. Ricco, V. K. Kozin, A. C. Seridonio, and I. A. Shelykh, Phys. Rev. A 106, 023702 (2022).
  * van Zanten  _et al._ [2020] D. M. T. van Zanten, D. Sabonis, J. Suter, J. I. Väyrynen, T. Karzig, D. I. Pikulin, E. C. T. O’Farrell, D. Razmadze, K. D. Petersson, P. Krogstrup, and C. M. Marcus, Nat. Phys. 16, 663 (2020).
  * Dai and Sun [2019] N. Dai and Q.-F. Sun, Phys. Rev. B 99, 085436 (2019).
  * Wang  _et al._ [2020] X.-Q. Wang, B. Wu, S.-F. Zhang, Q. Wang, and W.-J. Gong, Ann. Physics 415, 168127 (2020).
  * Wang  _et al._ [2021] X.-D. Wang, X.-Q. Wang, H.-Y. Sun, X.-F. Dai, Q. Wang, and W.-J. Gong, Phys. Scr. 96, 015805 (2021).
  * Wang [2021] Z.-H. Wang, Front. Phys. (Lausanne) 9, 704493 (2021).
  * Wang and Huang [2021] Z.-H. Wang and W.-C. Huang, Front. Phys. (Lausanne) 9, 727934 (2021).
  * Máthé  _et al._ [2022] L. Máthé, D. Sticlet, and L. P. Zârbo, Phys. Rev. B 105, 155409 (2022).
  * Zhu and Balatsky [2003] J.-X. Zhu and A. V. Balatsky, Phys. Rev. B 67, 165326 (2003).
  * Chen  _et al._ [2005] Z.-Z. Chen, R. Lü, and B.-f. Zhu, Phys. Rev. B 71, 165324 (2005).
  * Ramos-Andrade  _et al._ [2018] J. P. Ramos-Andrade, P. A. Orellana, and S. E. Ulloa, J. Phys.: Condens. Matter 30, 045301 (2018).
  * Jauho  _et al._ [1994] A.-P. Jauho, N. S. Wingreen, and Y. Meir, Phys. Rev. B 50, 5528 (1994).
  * Sun  _et al._ [2000] Q.-f. Sun, B.-g. Wang, J. Wang, and T.-h. Lin, Phys. Rev. B 61, 4754 (2000).
  * Michałek  _et al._ [2013] G. Michałek, B. R. Bułka, T. Domański, and K. I. Wysokiński, Phys. Rev. B 88, 155425 (2013).
  * Hwang  _et al._ [2015] S.-Y. Hwang, R. López, and D. Sánchez, Phys. Rev. B 91, 104518 (2015).
  * Bocian and Rudziński [2015] K. Bocian and W. Rudziński, Eur. Phys. J. B 88, 50 (2015).
  * Świrkowicz  _et al._ [2008] R. Świrkowicz, M. Wilczyński, and J. Barnaś, J. Phys.: Condens. Matter 20, 255219 (2008).
  * Máthé and Grosu [2018] L. Máthé and I. Grosu, Mater. Today:. Proc. 5, 15878 (2018).
  * Zhang  _et al._ [2012] S.-N. Zhang, W. Pei, T.-F. Fang, and Q.-f. Sun, Phys. Rev. B 86, 104513 (2012).
  * Máthé and Grosu [2020] L. Máthé and I. Grosu, Beilstein J Nanotechnol 11, 225 (2020).
  * Sun and Xie [2007] Q.-f. Sun and X. C. Xie, Phys. Rev. B 75, 155306 (2007).
  * Zimbovskaya [2008] N. A. Zimbovskaya, Phys. Rev. B 78, 035331 (2008).
  * Zhou  _et al._ [2019] L.-L. Zhou, X.-Y. Zhou, and R. Cheng, Chin. J. Phys. 60, 307 (2019).
  * Razmadze  _et al._ [2020] D. Razmadze, E. C. T. O’Farrell, P. Krogstrup, and C. M. Marcus, Phys. Rev. Lett. 125, 116803 (2020).


#  A First-Order Method with Expansive Projection for Optimal Powered Descent
Guidance

Jiwoo Choi  Graduate student, Department of Aerospace Engineering, Inha
University, Republic of Korea<EMAIL_ADDRESS>Jong-Han Kim  Associate
professor, Department of Aerospace Engineering, Inha University, Republic of
Korea<EMAIL_ADDRESS>(Corresponding author)

###### Abstract

This paper introduces a first-order method for solving optimal powered descent
guidance (PDG) problems, that directly handles the nonconvex constraints
associated with the maximum and minimum thrust bounds with varying mass and
the pointing angle constraints on thrust vectors. This issue has been
conventionally circumvented via lossless convexification (LCvx), which lifts a
nonconvex feasible set to a higher-dimensional convex set, and via linear
approximation of another nonconvex feasible set defined by exponential
functions. However, this approach sometimes results in an infeasible solution
when the solution obtained from the higher-dimensional space is projected back
to the original space, especially when the problem involves a nonoptimal time
of flight. Additionally, the Taylor series approximation introduces an
approximation error that grows with both flight time and deviation from the
reference trajectory. In this paper, we introduce a first-order approach that
makes use of orthogonal projections onto nonconvex sets, allowing expansive
projection (ExProj). We show that 1) this approach produces a feasible
solution with better performance even for the nonoptimal time of flight cases
for which conventional techniques fail and 2) the proposed method compensates
for the linearization error that arises from Taylor series approximation. We
claim that the proposed approach offers more flexibility in generating
feasible trajectories for a wide variety of planetary soft landing problems.

##  1 Introduction

Recently, the space industry has been undergoing a radical transformation with
the advent of the “new space” era. In particular, SpaceX’s Falcon 9 and Blue
Origin’s New Shepard have successfully demonstrated soft landing technologies
for reusable launch vehicles, and the European Space Agency is working on the
Themis program. In earlier days of planetary landing missions, classical
guidance algorithms such as those used in the Apollo guidance system were
developed, but they were limited to missions with various constraints.
Recently, convex optimization-based soft landing techniques have been
developed to handle a wide variety of practical constraints.

This paper addresses optimal powered descent guidance (PDG) problems with the
goal of finding minimum-fuel soft landing trajectories while considering a
variety of constraints on the state variables and control forces [1]. Two
notable constraints regarding PDG problems are as follows: 1) the magnitude of
the thrust is bounded not only by an upper limit but also by a lower limit
since the main engine cannot be turned off after ignition during the powered
descent phase, and 2) the thrust vector must lie within a certain range
because of maneuvering requirements or the field-of-view limits of the
navigation sensors. The lower bound constraints are obviously nonconvex, and
the pointing angle constraints can also be nonconvex depending on the range.
In previous works relying on techniques such as lossless convexification, the
problem has been reformulated by linearizing and relaxing these nonconvex
constraints by introducing a slack variable limiting the magnitude of the
thrust and the changes of variables handling mass variations [1, 2, 3].

It is known that the relaxation is exact, with the slack variable matching the
thrust magnitude, when the minimum fuel problem is solved for cases with the
optimal time of flight, $t_{f}^{*}$. However, when the problem is solved for a
nonoptimal time of flight, the solution obtained from the relaxation may be
infeasible for the original problem. Conventionally, the optimal time of
flight is found by applying bisection or similar search algorithms; however,
doing so increases the computational burden and hence is not desirable during
the powered descent phase, in which fast decision making is needed. Moreover,
the previous techniques additionally require Taylor series approximation for
handling the thrust bound constraints, which are reformulated in the form of
exponential functions. Such approximation tends to result in suboptimality or
constraint violation when the problem considers a significantly long time
horizon or when the obtained optimal solution deviates notably from the
reference trajectory along which the dynamics are linearized.

In this paper, we overcome the aforementioned difficulties by introducing a
first-order solution algorithm that directly handles nonconvex sets without
applying relaxation or linear approximation to the original PDG problem. It is
observed that for the _lossless convexification_ cases (with the optimal time
of flight), the proposed approach finds the same optimal solution that the
previously known approach finds. Additionally, the proposed approach is
observed to find very good practical solutions even for the _lossy
convexification_ cases (with a nonoptimal time of flight), for which the
previous approach fails.

We summarize the advantages of the proposed approach as follows:

  * •

It finds good practical solutions even when the previously known approach
fails.

  * •

It significantly reduces the suboptimality of the solution originating from
linear approximation error.

Notably, the proofs of convergence for first-order methods such as the
alternating direction method of multipliers (ADMM) or proximal gradient (PG)
techniques largely rely on the fact that projection onto convex sets is
_nonexpansive_. However, note that nonexpansivity of the projection operators
is not a necessary and sufficient condition for algorithm convergence; in
fact, it is an overly conservative condition, as there are many cases in which
expansive projections can still achieve global convergence. Our approach
involves orthogonal projection onto nonconvex sets, which can be _expansive_ ,
and is thus deserving of the name expansive projection (ExProj).

##  2 Powered Descent Guidance Problem

###  2.1 The PDG Problem in Nonconvex Form

The minimum-fuel PDG problem for planetary landers or reusable launchers is
defined as follows.

Problem (PDG):

$\displaystyle\underset{t_{f},T_{c}(\cdot)}{\text{minimize}}\quad$
$\displaystyle\displaystyle\int_{0}^{t_{f}}\left\|T(t)\right\|dt$ subject to
$\displaystyle\ddot{r}(t)=g+T_{c}(t)/m(t),$ (1a)
$\displaystyle\dot{m}(t)=-\alpha\left\|T(t)\right\|,$ (1b) $\displaystyle
0<\rho_{1}\leq\left\|T(t)\right\|\leq\rho_{2},$ (1c) $\displaystyle
e_{1}^{T}T(t)\geq\left\|T(t)\right\|\cos{\theta_{\text{tp}}},$ (1d)
$\displaystyle r(0)=r_{\text{init}},\ \dot{r}(0)=\dot{r}_{\text{init}},\
m(0)=m_{\text{wet}},$ (1e) $\displaystyle r(t_{f})=\dot{r}(t_{f})=0,\
m(t_{f})\geq m_{\text{dry}},$ (1f) $\displaystyle\forall t\in[0,t_{f}],$

where (1a) expresses the dynamics under the assumption that the vehicle is a
point mass and (1b) gives the change in mass due to thrust usage. The
constraints that follow, from (1c) to (1f), describe the limits on the thrust
magnitude, the pointing angle constraint, and the initial and final
constraints, respectively. In addition to the constraints described in Problem
(1), a wide range of additional constraints can be considered depending on the
mission and vehicle types. Note that the constraints in (1a) and (1b) are
nonlinear and that the lower bound on the thrust in (1c) is nonconvex. In
addition, (1d) can be nonconvex depending on $\theta_{\text{tp}}$. Moreover,
note that $t_{f}$ is an optimization variable rather than being fixed a
priori.

Table 1: Nomenclature $r(t)$ | Position vector
---|---
$\dot{r}(t)$ | Velocity vector
$T_{c}(t)$ | Thrust vector
$m(t)$ | Vehicle mass
$\rho_{1}$ | Minimum required thrust
$\rho_{2}$ | Maximum allowable thrust
$g$ | Gravitational acceleration vector
$\alpha$ | Constant describing mass consumption rate
$\theta_{\text{tp}}$ | Maximum thrust pointing angle from the vertical

###  2.2 Lossless Convexification

In the lossless convexification (LCvx) technique described in [1, 2, 3], a
slack variable $\Gamma(t)$ such that $\|T_{c}(t)\|\leq\Gamma(t)$ is introduced
to apply convex relaxation to the nonconvex constraints of Problem (1): the
lower bound on the thrust magnitude, $\rho_{1}\leq\left\|T_{c}(t)\right\|$,
and the pointing angle constraint on the thrust vector,
$e_{1}^{T}T_{c}(t)\geq\left\|T_{c}(t)\right\|\cos{\theta_{\text{tp}}}$.
Additionally, with the following changes of variables for linearizing the
dynamics,

$\displaystyle u(t)$ $\displaystyle=T(t)/m(t),$ $\displaystyle\sigma(t)$
$\displaystyle=\Gamma(t)/m(t),$ $\displaystyle z(t)$ $\displaystyle=\log
m(t),$ $\displaystyle z_{\text{ref}}(t)$
$\displaystyle=\log(m_{\text{wet}}-\alpha\rho_{2}t),$

the original constraint on the slack variable and the thrust bounds,
$\|T_{c}(t)\|\leq\Gamma(t)$ and $\rho_{1}\leq\Gamma(t)\leq\rho_{2}$,
respectively, are reformulated as follows.

$\|u(t)\|\leq\sigma(t),$ (2)
$\rho_{1}e^{-z(t)}\leq\sigma(t)\leq\rho_{2}e^{-z(t)}.$ (3)

Note that in (3), the inequality on the left-hand side defines a nonlinear but
convex set, while the inequality on the right-hand side defines a nonconvex
set. LCvx linearizes both sides around a reference trajectory
$z_{\text{ref}}(t)$, yielding the following convex problem.

Problem (LCvx):

$\displaystyle\underset{t_{f},\sigma(\cdot),u(\cdot)}{\text{minimize}}$
$\displaystyle\int_{0}^{t_{f}}\sigma(t)dt$ (4) subject to
$\displaystyle\ddot{r}(t)=g+u(t),$ $\displaystyle\dot{z}(t)=-\alpha\sigma(t),$
$\displaystyle\left\|u(t)\right\|\leq\sigma(t),$ $\displaystyle
e_{1}^{T}u(t)\geq\sigma(t)\cos{\theta}_{\text{tp}},$
$\displaystyle\rho_{1}e^{-z_{\text{ref}}(t)}\\{1-(z(t)-z_{\text{ref}}(t))\\}\leq\sigma(t),$
$\displaystyle\sigma(t)\leq\rho_{2}e^{-z_{\text{ref}}(t)}\\{1-(z(t)-z_{\text{ref}}(t))\\},$
$\displaystyle r(0)=r_{\text{init}},\ \dot{r}(0)=\dot{r}_{\text{init}},\
m(0)=m_{\text{wet}},$ $\displaystyle r(t_{f})=\dot{r}(t_{f})=0,\ m(t_{f})\geq
m_{\text{dry}},$ $\displaystyle\forall t\in[0,t_{f}].$

Note that Problem (4) is equivalent to Problem (1) only when the inequality in
(2) is tight.

$\|u(t)\|=\sigma(t).$ (5)

It has been proven that the optimal solution obtained from Problem (4) for the
optimal time of flight of Problem (1), $t_{f}^{*}$, is guaranteed to satisfy
$\left\|u^{*}(t)\right\|=\sigma^{*}(t),\forall t\in[0,t_{f}^{*}]$.

##  3 First-Order Optimization

###  3.1 First-Order Methods

First-order methods are a class of optimization algorithms that use the
gradient or subgradient information along with proximal operations to solve
optimization problems. Compared to second-order methods, which require the
Hessian information, first-order methods are computationally simpler and more
robust, easier to implement, and more efficient in terms of memory usage;
hence, they can handle problems on very large scales [4, 5, 11].

Although first-order methods can be less accurate and may require more
iterative computations than second-order methods, they are well suited for
control applications in which a high-accuracy solution is not necessarily
required and there are successive opportunities to update the solution in
subsequent time steps [9, 14].

###  3.2 Convergence, Nonexpansivity, and Expansive Projection

A general convex optimization problem with convex $f(\cdot)$ and $\mathcal{C}$
can be expressed as

$\underset{x\in\mathcal{C}}{\text{minimize}}\quad f(x)$

and can be solved via the ADMM, for example, by iteratively updating

$\displaystyle x^{k+1}$
$\displaystyle=\operatorname{prox}_{f/\rho}\left(z^{k}-u^{k}\right),$ (6)
$\displaystyle z^{k+1}$
$\displaystyle=\operatorname{prox}_{I_{\mathcal{C}}}\left(x^{k+1}+u^{k}\right),$
$\displaystyle u^{k+1}$ $\displaystyle=u^{k}+x^{k+1}-z^{k+1},$

where $I_{\mathcal{C}}(\cdot)$ denotes the indicator function of the set
$\mathcal{C}$.

The global convergence of (6) can be shown based on an understanding of the
nonexpansivity of the proximal operator of convex functions [4, 10, 11].
Noting that the orthogonal projection $\Pi_{\mathcal{C}}(\cdot)$ is the
proximal operator of the indicator function $I_{\mathcal{C}}(\cdot)$,
nonexpansivity for the convex $f(\cdot)$ and $\mathcal{C}$ implies that

$\begin{array}[]{c}\|\operatorname{prox}_{f}(x)-\operatorname{prox}_{f}(y)\|\leq\|x-y\|\quad\text{and}\vspace{1mm}\\\
\|\Pi_{\mathcal{C}}(x)-\Pi_{\mathcal{C}}(y)\|\leq\|x-y\|,\quad\forall
x,y.\end{array}$ (7)

On the other hand, as shown in Figure 1, the operator
$\Pi_{\mathcal{X}}(\cdot)$ for projection onto a nonconvex set $\mathcal{X}$
can be expansive.

$\left\|\Pi_{\mathcal{X}}(x)-\Pi_{\mathcal{X}}(y)\right\|\geq\left\|x-y\right\|,\quad\text{for
some }x,y.$

[width=0.6]plot

Figure 1: Orthogonal projections onto a nonconvex set $\mathcal{X}$. Observe
that projection onto nonconvex sets can be expansive.

Note that projection onto nonconvex sets can sometimes increase the distance
between two points, violating the nonexpansivity condition in (7). However,
this does not necessarily imply divergence of (6), as (7) is merely a
sufficient but not necessary condition for the convergence of (6). This
consideration can serve as a logical background for the intuitive application
of expansive projections to first-order optimization.

###  3.3 First-Order Method with Expansive Projection

The PDG problem convexified via LCvx can be readily solved using off-the-shelf
solvers. However, this approach is not able to produce a practically feasible
solution for some cases with $t_{f}$ not equal to $t_{f}^{*}$; moreover, the
solution may be suboptimal due to the approximation error arising from
linearizing (3).

In this paper, we do not apply convex relaxation or linear approximation.
Instead, we use a first-order method with direct projection onto nonconvex
sets, which can be expansive. With the nonconvex constraints in (3) and (5),
the problem can be defined as follows. Note that this problem is identical to
the original PDG problem in (1) and that two nonconvex constraints are present
in (8c) and (8e).

Problem (ExProj):

$\displaystyle\underset{t_{f},\sigma(\cdot),u(\cdot)}{\text{minimize}}\quad$
$\displaystyle-z(t_{f})$ subject to $\displaystyle\ddot{r}(t)=g+u(t),$ (8a)
$\displaystyle\dot{z}(t)=-\alpha\sigma(t),$ (8b)
$\displaystyle\left\|u(t)\right\|=\sigma(t),$ (8c) $\displaystyle
e_{1}^{T}u(t)\geq\sigma(t)\cos\theta_{\text{tp}},$ (8d)
$\displaystyle\rho_{1}e^{-z(t)}\leq\sigma(t)\leq\rho_{2}e^{-z(t)},$ (8e)
$\displaystyle r(0)=r_{\text{init}},\ \dot{r}(0)=\dot{r}_{\text{init}},\
m(0)=m_{\text{wet}},$ (8f) $\displaystyle r(t_{f})=\dot{r}(t_{f})=0,\
m(t_{f})\geq m_{\text{dry}}.$ (8g) $\displaystyle\forall t\in[0,t_{f}].$

###  3.4 ADMM with Expansive Projection

We apply the ADMM procedures with expansive projections (ExProj) onto the
nonconvex sets defined by (8c) and (8e). The ADMM algorithm combines dual
ascent and the method of multipliers to find the optimal solution by
alternately updating the primal variables and the dual variables [5]. With
$x=\begin{bmatrix}r^{T}&\dot{r}^{T}\end{bmatrix}^{T}$ and the terminal
constraints expressed in the soft constraint term, Problem (8) can be
discretized into the following standard form:

minimize $\displaystyle-z_{N}+\gamma\left\|x_{N}\right\|^{2}$ (9) subject to
$\displaystyle Py\geq q,$ $\displaystyle
y\in\mathcal{C}_{0}\cap\mathcal{C}_{1}\cap{C}_{2},$

with

$\displaystyle\mathcal{C}_{0}=\left\\{y\ |\ Gy=b\right\\},$ (10)
$\displaystyle\mathcal{C}_{1}=\cap_{i}\left\\{(u_{i},\sigma_{i})\ |\
\left\|u_{i}\right\|=\sigma_{i}\right\\},$
$\displaystyle\mathcal{C}_{2}=\cap_{i}\left\\{(z_{i},\sigma_{i})\ |\
\rho_{1}e^{-z_{i}}\leq\sigma_{i}\leq\rho_{2}e^{-z_{i}}\right\\},$

where $\gamma$ is some positive weighting parameter, and the horizon size $N$
satisfies $N\Delta t=t_{f}$ with sampling interval $\Delta t$. The vector
$y_{i}$ consists of the $i$-th state variables and control vector and takes
the form
$y_{i}=\begin{bmatrix}u_{i}^{T}&x_{i+1}^{T}&\sigma_{i}&z_{i+1}\end{bmatrix}^{T}$
for $i\in\left\\{0,1,\cdots,N-1\right\\}$, and we correspondingly define the
stacked variable
$y=\begin{bmatrix}y_{0}^{T}&y_{1}^{T}&\cdots&y_{N-1}^{T}\end{bmatrix}^{T}$.
Thus, we can encode the inequality constraints in (8d) as $Py\geq q$ and the
dynamic constraints in (8a) and (8b) as $Gy=b$. Accordingly, the problem is
further reformulated as follows:

minimize $\displaystyle-
z_{N}+\gamma\left\|x_{N}\right\|^{2}+I_{\mathcal{C}_{0}}(y)+I_{\mathcal{C}^{w}}(w)$
(11) subject to $\displaystyle D^{w_{1}}y-w_{1}=0,$ $\displaystyle
D^{w_{2}}y-w_{2}-b^{w_{2}}=0,$ $\displaystyle Py-q-w_{3}=0,$

with

$\displaystyle\mathcal{C}_{3}$ $\displaystyle=\left\\{w_{3}\ |\ w_{3}\geq
0\right\\},$ (12) $\displaystyle\mathcal{C}^{w}$ $\displaystyle=\left\\{w\ |\
w_{1}\in\mathcal{C}_{1},\ w_{2}\in\mathcal{C}_{2},\
w_{3}\in\mathcal{C}_{3}\right\\},$

where $w_{1}$ and $w_{2}$ are auxiliary variables of $y$ and $D^{w_{1}}$ and
$D^{w_{2}}$ are selection matrices that extract $w_{1}$ and $w_{2}$ from $y$.
The variable $w_{3}$ is a slack variable introduced for handling the
inequality constraints in $Py-q\geq 0$, and we additionally define the stacked
variable $w=\begin{bmatrix}w_{1}^{T}&w_{2}^{T}&w_{3}^{T}\end{bmatrix}^{T}$.

$\displaystyle D^{w_{1}}y$
$\displaystyle=\begin{bmatrix}u_{0}^{T}&\sigma_{0}&u_{1}^{T}&\sigma_{1}&\cdots&u_{N-1}^{T}&\sigma_{N-1}\end{bmatrix}^{T},$
$\displaystyle D^{w_{2}}y$
$\displaystyle=\begin{bmatrix}\sigma_{0}&\sigma_{1}&z_{1}&\sigma_{2}&z_{2}&\cdots&\sigma_{N-1}&z_{N-1}\end{bmatrix}^{T},$
$\displaystyle b^{w_{2}}y$
$\displaystyle=\begin{bmatrix}z_{0}&0&\cdots&0\end{bmatrix}^{T}.$

We define the augmented Lagrangian for Problem (11) as

$L_{\rho}(y,w,y_{s})=\overbrace{y^{T}Hy+hy+\frac{\rho}{2}\left\|Cy-\tilde{q}-w+y_{s}\right\|_{2}^{2}}^{f_{g}(y,w,y_{s})}+I_{\mathcal{C}_{0}}(y)+I_{\mathcal{C}^{w}}(w)$

with

$w=\begin{bmatrix}w_{1}\\\ w_{2}\\\
w_{3}\end{bmatrix},\,C=\begin{bmatrix}D^{w_{1}}\\\ D^{w_{2}}\\\
P\end{bmatrix},\,\tilde{q}=\begin{bmatrix}0\\\ b^{w_{2}}\\\
q\end{bmatrix},\,y_{s}=\begin{bmatrix}y_{s,1}\\\ y_{s,2}\\\
y_{s,3}\end{bmatrix},$

which can be solved by iteratively applying the following updates.

$\displaystyle y^{j+1}=$ $\displaystyle\
\underset{y}{\text{argmin}}\left(f_{g}(y,\ w^{j},\
y_{s}^{j})+I_{\mathcal{C}_{0}}(y)\right),$ $\displaystyle w^{j+1}=$
$\displaystyle\
\underset{w}{\text{argmin}}\left(\frac{\rho}{2}\left\|Cy^{j+1}-\tilde{q}-w+y_{s}^{j}\right\|^{2}+I_{\mathcal{C}^{w}}(w)\right),$
$\displaystyle y_{s}^{j+1}=$ $\displaystyle Cy^{j+1}-\tilde{q}-w^{j+1}.$

We handle the linear constraints $Gy=b$ without introducing an additional
multiplier. Instead, we make use of the KKT condition for the associated
equality-constrained quadratic problem [6].

$\displaystyle y^{j+1}=\ \underset{y}{\text{argmin}}\left(f_{g}(y,\ w^{j},\
y_{s}^{j})+I_{\mathcal{C}_{0}}(y)\right)\quad\Longleftrightarrow$
$\displaystyle\begin{bmatrix}\rho C^{T}C+H&G^{T}\\\ G&0\\\
\end{bmatrix}\begin{bmatrix}y^{j+1}\\\ *\end{bmatrix}=\begin{bmatrix}\rho
C^{T}\left(\tilde{q}+w^{j}-y_{s}^{j}\right)-h\vspace{5pt}\\\ b\end{bmatrix}.$

The updates for $w_{1}$ and $w_{2}$ involve a series of expansive projection
operations onto $\mathcal{C}_{1}$ and $\mathcal{C}_{2}$, respectively, while
the $w_{3}$ update is simple.

$\displaystyle w^{j+1}=$ $\displaystyle\
\underset{w}{\text{argmin}}\left(\frac{\rho}{2}\left\|Cy^{j+1}-\tilde{q}-w+y_{s}^{j}\right\|^{2}+I_{\mathcal{C}^{w}}(w)\right)$
$\displaystyle=$
$\displaystyle\begin{bmatrix}[l]\Pi_{\mathcal{C}_{1}}\left(D^{w_{1}}y^{j+1}-w_{1}+y^{j}_{s,1}\right)\\\
\Pi_{\mathcal{C}_{2}}\left(D^{w_{2}}y^{j+1}-w_{2}+y^{j}_{s,2}\right)\\\
\Pi_{\mathcal{C}_{3}}\left(Py^{j+1}-q+y^{j}_{s,3}\right)\end{bmatrix}.$

####  3.4.1 Projection onto $\mathcal{C}_{1}$

In general, projection onto the convex set defined by a second-order cone
constraint of the form $\left\|u\right\|\leq\sigma$ involves projection from
the exterior region of the cone [4, 14]. However, in our problem, projection
from both the exterior and interior regions should be considered because the
constraint in (8c) represents the surface of a second-order cone.

$\Pi_{\mathcal{C}_{1}}(u,\sigma)=\begin{cases}(u,\sigma)&\text{if
}\sigma=\|u\|\\\ 0&\text{if }\sigma\leq-\|u\|\\\
\left(\frac{\|u\|+\sigma}{2\|u\|}u,\frac{\|u\|+\sigma}{2}\right)&\text{otherwise}.\end{cases}$
(13)

Note from Figure 2 that the projection from region ③ can be expansive.

Figure 2: Orthogonal projection onto the _surface of the second-order cone_
defined by $\mathcal{C}_{1}$. Note that the set is nonconvex and that the
projection from region ① (exterior region) is nonexpansive, while the
projection from region ③ (interior region) can be expansive.

$1$$2$$3$$\sigma$

####  3.4.2 Projection onto $\mathcal{C}_{2}$

The set $\mathcal{C}_{2}$ describes an area surrounded by two exponential
curves. For a point $(z,\sigma)\notin\mathcal{C}_{2}$, the line passing
through the corresponding projected point and itself is perpendicular to the
slope at the projected point on the boundary. Based on this observation, we
can compute the projection as follows.

$\Pi_{\mathcal{C}_{2}}(z,\sigma)=\begin{cases}(z,\sigma)&\text{if
}(z,\sigma)\in\mathcal{C}_{2}\\\ (t_{1},\rho_{1}e^{-t_{1}})&\text{if
}\sigma<\rho_{1}e^{-z}\\\
(t_{2},\rho_{2}e^{-t_{2}})&\text{otherwise}.\end{cases}$ (14)

Here, $t_{1}$ and $t_{2}$ satisfy the followings.

$\displaystyle e^{t_{1}}(t_{1}-z)-\rho_{1}^{2}e^{-t_{1}}+\rho_{1}\sigma=0,$
(15) $\displaystyle
e^{t_{2}}(t_{2}-z)-\rho_{2}^{2}e^{-t_{2}}+\rho_{2}\sigma=0.$

Figure 3: Orthogonal projection onto $\mathcal{C}_{2}$ (shaded area). Note
that the set is nonconvex and that the projection from region ①
($\sigma<\rho_{1}e^{-z}$) is nonexpansive, while the projection from region ③
($\sigma>\rho_{2}e^{-z}$) can be expansive.

$1$$2$$3$

Note from Figure 3 that the projection from region ③ can be expansive.

Since the solution to (15) is not given in closed form, we numerically solve
the problem using the Newton–Raphson method, which rapidly converges within a
few steps for this problem; see Figure 4.

Figure 4: Relative residuals achieved with the Newton–Raphson method for the
instances from ExProj. Note that the calculation converges exponentially
within a few steps.

####  3.4.3 Projection onto $\mathcal{C}_{3}$

The projection onto the set $\mathcal{C}_{3}$ is simply obtained by taking
only the positive parts.

$\Pi_{\mathcal{C}_{3}}(Py^{j+1}-q+y^{j}_{s,3})=(Py^{j+1}-q+y^{j}_{s,3})_{+}.$

##  4 Numerical Examples

To verify the performance of the proposed approach (ExProj) and to make
quantitative comparisons with the previously known standard technique (LCvx),
we prepared three simulation cases. The first is the case of the optimal time
of flight, $t_{f}=t^{*}_{f}$, for which convexification with the LCvx
technique is lossless and both ExProj and LCvx are expected to find the
optimal solution. The other two cases concern nonoptimal time of flight,
$t_{f}<t^{*}_{f}$ and $t_{f}>t^{*}_{f}$, for which the convexification can be
lossy.

The simulation scenario is based on the Mars soft landing mission presented in
[1], and the simulation parameters are summarized in Table 2. The LCvx problem
was solved using the ECOS [8] solver via the cvxpy [7] parser, and the ExProj
problem was solved using a customized solver based on the algorithms presented
in this paper.

Table 2:  Simulation Parameters Parameter | Value |  | Parameter | Value
---|---|---|---|---
$r_{\text{init}}$ (m) | (2400, 450, -330) |  |  $m_{\text{wet}}$ (kg) | 2000
$\dot{r}_{\text{init}}$ (m/s) | (-10, -40, 10) |  |  $m_{\text{dry}}$ (kg) | 1700
$g$ (m/s2) | (-3.71, 0, 0) |  |  $\rho_{1}$ (N) | 4800
$\alpha$ (s/m) | 0.0005 |  |  $\rho_{2}$ (N) | 19200
$\theta_{\text{tp}}$ (deg) | 90 |  |  |

###  4.1 Optimal Time of Flight Case ($t_{f}=t^{*}_{f}$)

For the case with the optimal time of flight, $t_{f}=t^{*}_{f}$, the
convexification technique applied in LCvx is lossless, yielding
$\|u_{i}\|=\sigma_{i}$ for all $i=0,\dots,N-1$, and both LCvx and ExProj
converge to the globally optimal solution. The optimal time of flight for the
given scenario is found to be $t_{f}^{*}=46.96\ \text{s}$. The results are
summarized in Figure 5 and Table 3, which presents the final position and
velocity as well as the fuel consumption for the mission.

We observe that the two methods produce mostly the same results, with slight
differences that can be attributed to the approximation error due to the
linearization of (3) in LCvx. The effect of this approximation error is also
visible in the later part of the thrust profile, where the ExProj solution
utilizes the full maximum thrust, while the LCvx solution uses slightly less
than the maximum allowed thrust, making the LCvx solution slightly
conservative and suboptimal; see Table 3. However, note that the overall
effect of these differences on the achieved fuel consumption is not
significant in this case.

Figure 5:  Optimization results for the optimal time of flight case ($t_{f}=t^{*}_{f}$). Note that the convexification is lossless in this case, and both LCvx and ExProj find essentially the same optimal solution. However, note that in the later part of the thrust profile, the LCvx solution uses slightly less than the maximum allowed thrust. This makes the LCvx solution slightly suboptimal. Table 3: Optimization Results for $t_{f}=t_{f}^{*}$ | LCvx | ExProj
---|---|---
$r$ ($10^{-5}$m) | (0.978, -4.55, -3.93) | (0.959, -5.13, -3.89)
$\dot{r}$ ($10^{-1}$m/s) | (-0.367, -0.412, 1.67) | (-0.363, -0.355, 1.66)
$m_{\text{wet}}-m({t_{f}^{*}})$ (kg) | 201.00 | 200.66

###  4.2 Nonoptimal Time of Flight Cases ($t_{f}<t^{*}_{f}$ or
$t_{f}>t^{*}_{f}$)

Cases with $t_{f}\neq t_{f}^{*}$ are where the convexification in the LCvx
technique fails or the linearization error is significantly large.

For the first case, in which the time of flight is strictly shorter than the
optimal time ($t_{f}=41.8\,\text{sec}<t_{f}^{*}$), we observe that the
convexification due to LCvx fails with $\|u_{i}\|<\sigma_{i}$, resulting in an
infeasible solution. This can be seen in Figure 6, where the thrust profile
found with LCvx is below the lower limit, whereas the ExProj solution is still
feasible.

For the second case, in which the time of flight is strictly longer than the
optimal time ($t_{f}=82\,\text{sec}>t_{f}^{*}$), the convexification is
lossless with $\|u_{i}\|=\sigma_{i}$; however, the approximation error due to
the linear approximation of the thrust bound constraints in (3) means that the
LCvx approach considers incorrect limits, and the discrepancy increases as the
flight time increases. This is clearly observed from the last part of the
thrust profile in Figure 7, where the LCvx solution uses significantly less
thrust than the allowed maximum. Consequently, the LCvx solution is
significantly conservative and thus suboptimal, whereas the ExProj solution is
free of this suboptimality.

Figure 6: Optimization results for the case in which the time of flight is
strictly shorter than the optimal time ($t_{f}<t_{f}^{*}$). Note that in this
case, the convexification technique from the LCvx approach is no longer
lossless, and the obtained solution violates the lower bound limit on the
thrust. On the other hand, ExProj successfully finds a feasible solution.

Figure 7: Optimization results for the case in which the time of flight is
strictly longer than the optimal time ($t_{f}>t_{f}^{*}$). Note that the
effect of the linear approximation error on the thrust profile increases with
the flight time. This is due to the linear approximation of (3) used in the
LCvx technique, which makes the LCvx solution suboptimal.

##  5 Experimental Validation

We briefly present the results obtained from flight tests conducted using a
drone in the indoor flight arena equipped with motion capture systems. The
proposed algorithm was implemented on a GPU to efficiently compute and update
the optimal PDG trajectory at a rate of 10 computations per second.
Subsequently, the computed thrust command is transformed into both throttle
and attitude commands of the drone. In this test, we configured the parameter
$\alpha$ in (1b) to have a very small value, as the weight of the drone
remains constant regardless of thrust usage. Furthermore, we incorporated a
dynamic adjustment to the maximum allowable tilt angle of the thrust vector,
which pregressively decreases as the vehicle approaches the landing pad [12].

Figure 8 presents the trajectory and thrust vector data acquired from the
flight test, along with the corresponding tilt angle limit. The figure also
includes the LCvx solution (depicted by lighter lines) computed from the
identical initial conditions. Notably, it is observed that the LCvx solution
is infeasible, necessitating a thrust level below the estabilished lower
bound. Finally, in Figure 9, we present the trajectory and a sequence shot
obtained through the motion capture systems.

(a) Trajectory.

StartFinish

(b) Thrust and the tilt angle.

Figure 8: Trajectory and the thrust profile from the indoor flight test. The
flight results obtained through ExProj (depicted by solid lines) are compared
to the computational results produced by LCvx (depicted by lighter lines),
computed using the identical conditions. Note that the thrust profile derived
from ExProj is feasible.

(a) Trajectory obtained from the motion capture system.

(b) Sequence shot from the flight test.

StartFinish

Figure 9: Experimental validation of ExProj via indoor flight test.

##  6 Concluding Remarks

In this paper, we proposed a first-order method that directly handles the
nonconvexity arising in powered descent guidance problems. Our approach
combines first-order convex optimization algorithms with orthogonal
projections onto nonconvex sets, which are sometimes expansive. Through a
series of numerical examples, we verified the performance of the proposed
algorithm and compared it with the most well-known standard convexification
approach.

The results of the numerical experiments show that our approach computes the
same optimal solution as the standard convexification technique when the
latter is lossless. More importantly, our approach finds practical feasible
solutions even when the standard technique fails and yields superior solutions
when the solutions obtained with the standard technique are significantly
suboptimal. Also, we briefly presented results from an indoor flight test
demonstrating the effectiveness of the proposed algorithm.

Although some results on the convergence of first-order methods on nonconvex
problems have been reported recently [13], general results are not yet known.
Therefore, a natural extension of the present work will be to analyze the
convergence of this specific algorithm.

## References

  * [1] B. Açıkmeşe, J. M. Carson, and L. Blackmore.  Lossless convexification of nonconvex control bound and pointing constraints of the soft landing optimal control problem.  IEEE Transactions on Control Systems Technology, 21(6):2104–2113, 2013.
  * [2] B. Açıkmeşe and S. R. Ploen.  Convex programming approach to powered descent guidance for mars landing.  Journal of Guidance, Control, and Dynamics, 30(5):1353–1366, 2007\.
  * [3] B. Açıkmeşe, D. Scharf, L. Blackmore, and A. Wolf.  Enhancements on the convex programming based powered descent guidance algorithm for mars landing.  In AIAA/AAS astrodynamics specialist conference and exhibit, page 6426, 2008.
  * [4] A. Beck.  First-order methods in optimization.  SIAM, 2017.
  * [5] S. Boyd, N. Parikh, E. Chu, B. Peleato, J. Eckstein, et al.  Distributed optimization and statistical learning via the alternating direction method of multipliers.  Foundations and Trends® in Machine learning, 3(1):1–122, 2011.
  * [6] M. S. Darup.  Encrypted mpc based on ADMM real-time iterations.  IFAC-PapersOnLine, 53(2):3508–3514, 2020.
  * [7] S. Diamond and S. Boyd.  CVXPY: A Python-embedded modeling language for convex optimization.  Journal of Machine Learning Research, 17(83):1–5, 2016.
  * [8] A. Domahidi, E. Chu, and S. Boyd.  ECOS: An SOCP solver for embedded systems.  In European Control Conference (ECC), pages 3071–3076, 2013.
  * [9] P. Elango, A. G. Kamath, Y. Yu, B. Açıkmeşe, M. Mesbahi, and J. M. Carson.  A customized first-order solver for real-time powered-descent guidance.  In AIAA SciTech 2022 Forum, page 0951, 2022.
  * [10] N. Parikh, S. Boyd, et al.  Proximal algorithms.  Foundations and trends® in Optimization, 1(3):127–239, 2014.
  * [11] E. K. Ryu and W. Yin.  Large-scale convex optimization: algorithms & analyses via monotone operators.  Cambridge University Press, 2022.
  * [12] D. P. Scharf, M. W. Regehr, G. M. Vaughan, J. Benito, H. Ansari, M. Aung, A. Johnson, J. Casoliva, S. Mohan, D. Dueri, et al.  ADAPT demonstrations of onboard large-divert guidance with a vtvl rocket.  In 2014 IEEE aerospace conference, pages 1–18. IEEE, 2014.
  * [13] Y. Wang, W. Yin, and J. Zeng.  Global convergence of admm in nonconvex nonsmooth optimization.  Journal of Scientific Computing, 78:29–63, 2019.
  * [14] Y. Yu, P. Elango, and B. Açıkmeşe.  Proportional-integral projected gradient method for model predictive control.  IEEE Control Systems Letters, 5(6):2174–2179, 2020.