\section{\label{sec:level1}introduction}
	Transition metal oxides (TMOs) exhibit diverse physical phenomena such as  metal-insulator transition \cite{ImadaFujimori_MIT_RMP1998}, superconductivity \cite{Bednorz1986}, multiferroicity \cite{M.Dawber_Ferroelectricoxides_RevModPhys.77.1083}, giant and colossal magnetoresistance \cite{ M.B.Salamon_GMR_RevModPhys.2001}, non-Fermi liquid behaviour \cite{Pkhalifah, PKhalifah_CSRO_NFLPhysRevB2004}, quantum phase transition \cite{Z.Q.Mao_BaRu6O12_PRL2003} and various exotic magnetic orders \cite{PBAllen_SrRuO3_PhysRevB1996, RSPErry_MetamagnetismSr3Ru2O7_PRL2001, S.A.Grigera_Sr3Ru2O7_Science2004,C.-Q.Jin_FerromagnetismARuO3_PNAS2008, S.Trebst_Kitaev_PhysicsReports20221}. It is well recognized that electron correlation plays a crucial role in describing such exotic properties. The electron correlation is expected to be weak in 4$d$ TMOs compared to 3$d$ TMOs due to larger spatial extension of 4$d$ orbitals than that of 3$d$ orbitals. However, against this general belief, varying strength of electron correlation have been observed in 4$ d $ TMOs leading to exotic ground states \cite{G.Koster_SRO_RevModPhys2012, HDKim_Ru3d_PhysRevLett2004, S.Nakatsuji_CaSrRuO4_PhysRevLett.2000, Maeno1994, Y.Park_SrNbO3_CommPhys2020,B.Chiara_SrNbO_ARPES_PRM2020, KMaiti_CaSaRuO3_PRB2005,KarpJSr2MoO42020, YZhangSrRhO3PRB2020,M.W.Havekort_Sr2RhO4._PhysRevLett2008, BHReddy_APd3O4_JPCM2021}.

	Among 4$d$ TMOs, molybdenum based perovskites $A$MoO$_{3}$ ($A$ = Ca, Sr, Ba) exhibit metallic behaviour with Pauli paramagnetism \cite{HAYASHI_AMoO3_1979}. Exceptionally high conductivity is found in SrMoO$_{3}$ associated with highly delocalized Mo 4$d^{2}$ electrons with weak contribution from phonon scattering \cite{NagaiAPL2005}. Resistivity measurement shows a Fermi-liquid behaviour and the Sommerfeld coefficient evaluated from specific heat measurements is twice than that obtained from band structure calculation, suggesting an enhanced quasiparticle mass \cite{NagaiAPL2005}. Such quasiparticle mass enhancement is considered as signature of electron correlation, also observed in other strongly correlated systems \cite{Inoue_CSVO_PhysRevB1998, PBAllen_SrRuO3_PhysRevB1996}. SrMoO$_{3} $ exhibits two structural transitions with lowering temperature where the room temperature cubic $Pm\bar{3}m$ structure goes to tetragonal $I4/mcm$ structure and further to orthorhombic $Imma$ structure at around 266 K and 124 K, respectively \cite{Macquart2010}. The low temperature pseudo-cubic structures arise due to the rotation and/or tilting of MoO$_{6}$ octahedra. Phonon calculations within DFT + $U$ ($U$: on-site Coulomb repulsion) framework emphasizes that the electron correlation is responsible for the structural transitions \cite{J.-L.Hand_SrMoO3_PRM2021}. Also, DFT + DMFT correctly predicts the non-magnetic ground state along with octahedra rotation for the $ Imma $ structure \cite{AHampel_SrMoO3_PRB2021}. Quasiparticle spectral weight ($ Z $) is found to be $\sim$ 0.6 \cite{AHampel_SrMoO3_PRB2021, APaulPRB2019}, commensurate with specific heat measurements and recent angle-resolved photoemission spectroscopy (ARPES) measurement \cite{NagaiAPL2005, E.Cappelli_SrMoO3ARPES_PRM2022}. However, DFT + DMFT fails to capture an intense satellite observed at $\sim$ 2.5 eV binding energy in hard $ x $-ray valence band photoemission spectra. This satellite has been argued to have a plasmonic origin \cite{HWadatiSrMoO3}. Various GW calculations have also predicted a plasmon satellite at about 3 eV binding energy but with much smaller spectral weight in contrast to the experiment \cite{NilssonPRM2017, PetochhiPRR2020, TianyuCellDMFTPRX2021}.

	In this letter, we investigate the electronic structure of SrMoO$_{3}$ using valence band photoemission spectroscopy on \textit{in-situ} fractured polycrystalline sample. We observe considerably different surface and bulk electronic structure. Observation of significantly \textit{weak} plasmon satellite along with moderate electron correlation in the bulk Mo 4$d$ spectra, are consistent with many-body theoretical calculations. We further discuss the surface electronic structure.

	High quality polycrystalline sample of SrMoO${_3}$ was prepared by solid state reaction method using high purity MoO$_{3}$ (99.99 \%) and SrCO$_{3}$ (99.995 \%). Thoroughly ground mixture was palletized and heated at 600 $^{\circ}$C for 8 hours and further at 1250 $^{\circ}$C for 48 hours with intermittent grindings. Heat treatments were performed under 5$\%$ hydrogen mixed argon gas flow, resulting in very hard brick-red pellets. The phase purity and crystal structure were confirmed by \textit{x}-ray diffraction pattern collected at room temperature (Fig. S1 of supplemental material (SM) \cite{Supplemental}). The cubic lattice parameter was found to be 3.975(8) \AA ~in excellent agreement with earlier reports \cite{NagaiAPL2005,Macquart2010}. Photoemission spectroscopic measurements were carried out at 30 K using monochromatic Al $K_{\alpha}$ ($h\nu$ = 1486.6 eV) and He {\scriptsize II} ($h\nu$ = 40.8 eV) radiation (energy) on \textit{in-situ} fractured sample (base pressure $\sim$ 4$\times$10$^{-11}$ mbar). The Fermi level ($E_{F}$) and energy resolution were determined by measuring the Fermi cut-off of a clean polycrystalline silver at 30 K. The energy resolutions for Al $K_{\alpha}$ and He {\scriptsize II} spectra were set to $\sim$ 300 meV and $\sim$ 10 meV, respectively.

	Electronic structure calculations were performed for the orthorhombic \textit{Imma} structure with structural parameters adopted from Ref. \cite{Macquart2010}. Full potential linearized augmented plane wave method, as implemented in \textsc{wien2k} \cite{WIEN2k}, was used for the DFT calculation within generalized gradient approximation of Purdew-Burke-Ernzerhof \cite{PBEPRL1996}. For the DFT-hybrid calculations, screened hybrid functional (YS-PBE0) \cite{TranYSPBE0PRB2011} was constructed by replacing the $\alpha$ fraction of semi-local exchange with Hartree-Fock exchange, where $\alpha$ is the mixing parameter. The energy and charge convergence criteria were set to $10^{-5}$ eV and $10^{-4}$ electronic charge per formula unit (f.u.), respectively.
	Fully charge self-consistent DFT + DMFT calculation for 50 K ($\beta$ $\approx$ 232 eV$^{-1}$) were performed using eDMFT code with continuous time quantum Monte Carlo impurity solver and ``exact" double counting \cite{HauleDMFTPRB2010,Haule_CTQMC_PRB2007, HauleExactDCPRL2015}. A hybridization window of $\pm$ 10 eV was used and all five $d$ orbitals as basis was chosen for correlated Mo atoms. Analytical continuation was performed using maximum entropy method to calculate self-energy on the real axis \cite{HauleDMFTPRB2010}. The chosen local axes were nearly aligned with Mo-O bond direction of the MoO$_{6}$ octahedra. The Hubbard-$U$ and Hund's coupling $J$ were set to 6.0 eV and 0.7 eV, respectively. Total 4000 $k$-points were used for DFT and DFT + DMFT calculations and 500 $k$-points were used for DFT-hybrid calculation. Density of states (DOS) was calculated using 4000 $k$-points for all the calculations.

\begin{figure}[t]
	\centerline{\includegraphics[width=.45\textwidth]{Fig1.eps}}
	\caption{(color online) (a) Valence band photoemission spectra of SrMoO$_{3}$ collected at 30 K using Al $K_{\alpha}$ (black circles) and He \textsc{ii} (red circles) radiations. The O 2$p$ band contributions (see text) are shown by blue line. Atomic arrangement in orthorhombic structure is shown where red, black, pink and blue spheres represents Sr, Mo, O1 (apical) and O2 (basal) atoms, respectively. (b) DFT calculated PDOS for Sr, Mo, O1 and O2 atoms for the orthorhombic $Imma$ structure.}\label{fig:Fig1}
\end{figure}

	Atomic arrangement in orthorhombic SrMoO$ _{3} $ has been shown in Fig \ref{fig:Fig1} (a). Each Mo atom is surrounded by two types of  crystallographically different oxygens (two apical (O1) and four basal (O2)), while the Sr atom sits in the void created by surrounding MoO$ _{6} $ octahedra. Valence band is formed by the hybridization between Mo 4$d$ and O 2$p$ states and contributions from Sr states are expected to be negligible in the occupied energy range. Valence band photoemission spectra of SrMoO$_{3}$ collected at 30 K using Al $K_{\alpha}$ and He \textsc{ii} radiations have been shown in Fig. \ref{fig:Fig1} (a). Both spectra exhibit two distinctly separated groups of features below and above 3.5 eV binding energy. Considering the larger photo-ionization cross-section ratio of Mo 4$d$ states with O 2$p$ states in case of Al $K_{\alpha}$ than in case of He \textsc{ii} \cite{YEH19851}, the features below 3.5 eV binding energy can be attributed to Mo 4$d$ states and the features above 3.5 eV binding energy can be attributed to O 2$p$ states.

	The attributed characters of the spectral features are further examined using DFT calculation. Calculated partial density of states (PDOS) for Sr, Mo, O1 and O2 atoms are shown in Fig. \ref{fig:Fig1} (b). It is clear that the valence band is formed by Mo 4$d$ and O 2$p$ hybridized states, while Sr states have negligible contribution in the occupied region (appearing between -2 eV to -8 eV). There are three sets of features above $\sim$ 3.5 eV binding energy centered $\sim$ 4.5 eV corresponding to non-bonding states primarily having O 2$ p $ character and $\sim$ 6 eV and $\sim$ 8 eV features corresponding to bonding states with Mo 4$ d $ and O 2$ p $ mixed character. The anti-bonding states located between 2 eV to -2 eV binding energy primarily have Mo 4$d$ $t_{2g}$ character, while the Mo 4$d$ $e_{g}$ states appear between -2 eV to -7.2 eV binding energy. The average $\angle$ Mo-O-Mo reduces to $\sim$ 172.6$^{\circ}$ ($\angle$Mo-O1-Mo = $\sim$ 171.3$^{\circ}$ and $\angle$Mo-O2-Mo = $\sim$ 173.2$^{\circ}$) in the orthorhombic $Imma$ structure from 180$^{\circ}$ in the cubic $ Pm\bar{3}m$ structure \cite{Macquart2010}. It is to note here that the PDOS corresponding to O1 and O2 ($\times$ 0.5) are almost degenerate, as evident from the figure, suggesting negligible influence of the octahedra rotation and/or tilt in the electronic structure of SrMoO$_{3}$. This is also confirmed by very similar calculated DOS and essentially similar valence band spectra collected at 300 K and 30 K using Al $K_{\alpha}$, (Fig. S2 of SM \cite{Supplemental}).

\begin{figure}[t]
	\centerline{\includegraphics[width=.45\textwidth]{Fig2.eps}}
	\caption{(color online) Total DOS of orthorhombic SrMoO$_{3}$ calculated using DFT, DFT-hybrid and DFT + DMFT.}\label{fig:Fig2}
\end{figure}

	The overall comparison of experimental spectra with DFT results reveal two significant differences, (i) a mismatch in the energy position of the O 2$p$ band and (ii) an additional satellite feature around 3 eV binding energy in the experimental spectra. The feature at $\sim$ 5 eV binding energy in the experimental spectra appears at about 4.5 eV in the DOS calculated within DFT. This overestimation of O 2$p$ band energies has been attributed to the underestimation of electron correlation in DFT, as also observed in various TMOs \cite{SarmaDD_PRL1995}.
	The methods beyond DFT, such as DFT-hybrid and DFT + DMFT, have been found to be quite successful in case of TMOs with varying strength of electron correlation. The hybrid functionals containing some part of the Hartree-Fock exchange have been very successful in correctly describing the electronic properties of many semiconductors \cite{TranYSPBE0PRB2011,HSE2004, BHReddy_APd3O4_JPCM2021}, insulators \cite{TranYSPBE0PRB2011, JPaier_hybriddielectric_PRBR2008, SBansal_Li2RuO3_JPCM2021} and correlated metallic systems \cite{MCasedei_Cerium_PRB2016, Safdat_PbMoO3_2017, FalkePRB2021}. The mixing parameter $\alpha$ is related to the dielectric properties of a material and have been varied in order to match with the experimental results \cite{Koller_JPCM2013} even in case of metals \cite{FalkePRB2021}. Indeed, a shift of about 0.5 eV towards higher binding energy in the O 2$p$ band position is found for $\alpha$ = 0.10 as shown in the Fig. \ref{fig:Fig2}. DFT-hybrid calculations for cubic structure with varying $\alpha$ from 0.25 to 0.10 has been shown in SM \cite{Supplemental}. To further examine the observed shift, we have also performed DFT + DMFT calculation. The DFT + DMFT has emerged as a successful method to treat weak to strong electron correlation in many $d$ and $f$ electron systems \cite{GKotliar_RevModPhys.78.865, G.Kotliar_PhysicsToday2004}. The total DOS calculated using DFT + DMFT has been shown in Fig. \ref{fig:Fig2} which also exhibits shift of the O 2$p$ band towards higher binding energy as compared to DFT, consistent with the experimental spectra and DFT-hybrid. Reduced Mo 4$d$ bandwidth in DFT + DMFT suggests strong renormalization with quasiparticle weight $Z$ $\approx$ 0.5 (for $t_{2g}$ orbitals), consistent with earlier calculations and specific heat measurements \cite{APaulPRB2019, AHampel_SrMoO3_PRB2021, E.Cappelli_SrMoO3ARPES_PRM2022, NagaiAPL2005}. The DFT + DMFT provides a better description of the experimental spectra, as discussed later.

\begin{figure}[t]
	\centerline{\includegraphics[width=.45\textwidth]{Fig3.eps}}
	\caption{(color online) Extracted Mo 4$d$ band of SrMoO$ _{3} $ from Al $K_{\alpha}$ (black circles) and He \textsc{ii} (red circles) valence band spectra. The resolution broadened He \textsc{ii} spectra is shown by red line.}\label{fig:Fig3}
\end{figure}

	Interestingly, all of the calculations discussed above fails to describe the $\sim$ 3 eV satellite feature of the experimental spectra. Such higher binding energy feature observed in the photoemission spectra is a typical signature of correlation induced lower Hubbard band (LHB) in strongly correlated TMOs and has been successfully captured by DFT + DMFT \cite{A.Sekiyama_CaSrVO_PRL2004,Silke_VO2cDMFT_PRL2005, Z.V.Pchelkina_LRODMFT_PRB2015}. However, in the present case the satellite feature can not be attributed to LHB since it remains absent in the DFT + DMFT. A similar feature observed in earlier photoemission experiment was suggested to be a plasmon satellite \cite{HWadatiSrMoO3}, as also observed in subsequent GW calculations \cite{NilssonPRM2017, PetochhiPRR2020,TianyuCellDMFTPRX2021}.

\begin{figure}[t]
	\centerline{\includegraphics[width=.45\textwidth]{Fig4.eps}}
	\caption{(color online) Upper panel shows bulk spectra (open circles) of SrMoO$_{3}$. The solid blue and dashed pink lines shows the DFT + DMFT spectral function and Gaussian type bulk plasmon satellite (P$ _{B} $), respectively. Lower panel shows surface spectra (open circles) of SrMoO$ _{3} $. The solid grey, dashed green and dashed maroon lines shows bulk coherent feature, surface plasmon satellite (P$ _{S} $) and surface LHB, respectively.}\label{fig:fig4}
\end{figure}

	In order to further investigate, we extract the Mo 4$d$ band contributions from the valence band spectra. This can be reliably done by subtracting the O 2$p$ band features (simulated using three Gaussian and shown as blue lines), since they are distinctly separate from Mo 4$d$ band in experimental and calculated valence band shown in Fig. \ref{fig:Fig1}. The extracted Mo 4$d$ bands, normalized by total integrated intensity, are shown in Fig. \ref{fig:Fig3}, for both Al $K_{\alpha}$ and He \textsc{ii} spectra. Both spectra exhibit intense feature close to $E_{F}$ with a hump like satellite feature. The relatively weak satellite feature observed here is consistent with our previous $x$-ray photoemission study on \textit{ex-situ} thin film \cite{AAli_SMOfilm_AIPConf2019}. He \textsc{ii} spectra exhibits a sharper Fermi cut-off due to higher resolution than that of Al $K_{\alpha}$ spectra. For a direct comparison, we broaden the He \textsc{ii} spectra by a Gaussian of 0.3 eV width ($\sim$ energy resolution of Al $K_{\alpha} $ spectra) as shown by line. Now, these two spectra having similar resolution broadening differ from each other only in probing depth employed in the photoemission spectroscopy where Al $K_{\alpha}$ spectra is more bulk sensitive while He \textsc{ii} spectra is more surface sensitive. Different line shape of these two spectra suggests that the surface and bulk electronic structures are distinctly different in this system. Lowered symmetry, reduced coordination number and/or surface reconstruction $ etc.$ at the surface may lead to difference in the electronic structure. Thus, it is essential to disentangle the contribution of the surface to understand the intrinsic bulk electronic structure \cite{KMaiti_LACaVO_PRL1998,KMaiti_CaSaRuO3_PRB2005}.

	Photoemission spectral intensity for the incident photon energy $h\nu$ can be expressed as $I_{h\nu}(E) = e^{-l/\lambda_{h\nu}}f_{b}(E)+ (1-e^{-l/\lambda_{h\nu}})f_{s}(E)$, where $l$ is surface layer thickness, $\lambda$ is photoelectron mean free path (probing depth) and $f_{b}(E)$ and $f_{s}(E)$ represent bulk and surface spectra, respectively. $f_{b}(E)$ and $f_{s}(E)$ were estimated using $l/\lambda_{Al K_{\alpha}}$ = 0.45 and $l/\lambda_{He \textsc{ii}}$ = 1.75 and have been shown in Fig. \ref{fig:fig4}. The obtained spectra are quite robust with respect to the energy positions and relative intensity of features within 20$\%$ variation of $l/\lambda$, providing confidence in the analysis \cite{KMaiti_LACaVO_PRL1998,KMaiti_CaSaRuO3_PRB2005}.

	As evident from Fig. \ref{fig:fig4}, surface and bulk spectra (normalized by total integrated intensity) are distinctly different. The intrinsic bulk spectra exhibiting an intense coherent feature below 2 eV binding energy is compared with resolution broadened occupied part of the DFT + DMFT spectral function (matched at highest intensity). An excellent agreement for the coherent feature confirms moderate electron correlation in this system. A significantly \textit{weak} satellite (P$ _{B} $) in the bulk spectra at $\sim$ 3 eV binding energy is commensurate with the broad and weak plasmon satellite observed in various GW calculations \cite{NilssonPRM2017, PetochhiPRR2020, TianyuCellDMFTPRX2021}. It is to note here that the DFT, DFT-hybrid and DFT + DMFT calculations, performed here, fail to describe the plasmonic satellite, since, these calculations do not incorporate any long-range (non-local) interactions which are essential for description of plasmonic excitations \cite{HWadatiSrMoO3,Boehnke_PRB2016}.

	The surface spectra exhibits significantly reduced total width and enhanced satellite feature (appearing below 2.5 eV) in contrast to the bulk spectra. Interestingly, the surface spectra obtained here is strikingly similar to hard $x$-ray photoemission spectra on \textit{ex-situ} thin films \cite{HWadatiSrMoO3}. Plasmon satellites in the core level photoemission have been extensively studied revealing different lineshape and intensity for the bulk and surface plasmon. Surface plasmon satellite is expected to appear at smaller energy ( $\sqrt{2}$ times lower) than bulk plasmon satellite due to reduced dimensionality and/or confinement effects \cite{R.H.Ritchie_PhysRev.106.874_1956, Inglesfield_1983}. These have indeed been observed in grazing angle core-level photoemission experiments \cite{Biswas_AlPlasmon_PRB2003} while the enhancement of surface plasmon satellites are attributed to the enhancement of extrinsic plasmon losses at the surface \cite{Yubero_calculationAlplasmon_PRB2005, YUBERO_Siplasmon_SUrfScience2005}. A careful look at the surface spectra reveals that the coherent feature remains at very similar position while an intense and shifted satellite suggests an enhanced surface plasmon (P$ _{S} $), roughly $\sqrt{2}$ times lower energy than the bulk plasmon (P$ _{B} $). Deconvolution of the surface spectra as shown in Fig. \ref{fig:fig4}, requires at least two Gaussian type features in addition to bulk coherent feature (obtained from bulk spectra after subtracting Gaussian type P$ _{B} $ and matched at highest intensity), suggesting an additional feature at $\sim$ 1.5 eV binding energy (Gaussian peaks have been used for simplistic illustration). This additional feature can be a signature of LHB appearing due to enhanced effective electron correlation at the surface and requires further investigations on a high-quality crystals/thin films.

	In conclusion, we have investigated the electronic structure of SrMoO$_{3}$ using valence band photoemission spectroscopy. An accurate description of the O 2$p$ band position is found with DFT-hybrid and DFT + DMFT calculations. Valence band obtained using different photon sources suggest that the surface and bulk electronic structures are quite different in this system. The large coherent feature in bulk spectra is commensurate with DFT + DMFT spectral function suggesting moderate electron correlation in the Mo 4$d$ orbitals. Enhanced plasmon satellite and signature of strong correlation induced LHB are observed in the surface spectra. Intrinsic bulk spectra exhibits a \textit{weak} plasmon satellite as also observed in various GW calculations.

	We thank S. K. Pandey (IIT Mandi) for fruitful discussions. We acknowledge the support of CIF and HPC facilities at IISER Bhopal. We also thankfully acknowledge the funding from DST-FIST (Project No. SR/FST/PSI-195/2014C).


\section{Introduction}\label{sec1}

Let $(Z,d)$ be a metric space and $CL(Z)$ denotes the collection of non-void closed subsets of $Z$. The set $CB(Z)$ contains all the non-void bounded closed subsets of $Z$. For any $P,Q\in CB(Z)$,
\begin{equation}
	H(P,Q)=\max\bigg\{\sup_{u\in Q}D(u,P), \sup_{v\in P}D(v,Q)\bigg\}
\end{equation}
is called Hausdorff metric on $CB(Z)$, where $D(u,P)=\inf\{d(u,v):v\in P\}$. For a mapping $F:Z\to CL(Z)$, an element $z^*\in Z$ is known as a fixed point for the map $F$ if $z^*\in Fz^*$. The occurrence of fixed points for set-valued mappings was first deduced by Nadler \cite{nadler} in the year of 1969. Later, Mizoguchi and Takahashi \cite{mizo} improved the Nadler's \cite{nadler} theorem for a set-valued generalized contraction in the year 1989. The authors \cite{mizo} established the below stated theorem.
\begin{theorem}\cite{mizo}\label{mizo_thm}
	Suppose that $(Z,d)$ is complete and a set-valued map $F:Z\to CB(Z)$ in order that for every $y,z\in Z$,
	\begin{equation}\label{mizo_eq}
		H(Fy,Fz)\leq k(d(y,z))d(y,z),
	\end{equation}
	where $k\in W=\{\theta:[0,\infty)\to [0,1)| \limsup_{t\to s^+}\theta(t)<1,~\forall s\in [0,\infty)\}$. Then there is an element $z^*\in Fz^*$.
\end{theorem}
In fact, Suzuki \cite{suzuki} provided an example to deduce that every Mizoguchi -Takahashi contraction \cite{mizo} is not necessarily a Nadler contraction \cite{nadler} on arbitrary metric spaces in the year of 2007. However, Eldred-Anuradha-Veeramani \cite{eldred} observed that any Mizoguchi-Takahashi \cite{mizo} contraction follows a Nadler \cite{nadler} contraction, whenever the metric space $(Z,d)$ is metrically convex complete \cite{eldred} in the year 2009. There are several extensions of the Theorem \ref{mizo_thm} due to Mizoguchi and Takahashi \cite{mizo} in the literature, which can be found in \cite{berinde,mam}.

In the year 1961, Edelstein \cite{edel} generalized the famous Banach contraction principle for uniformly local contraction \cite{edel}. Subsequently, Nadler \cite{nadler} improved the Edelstein's \cite{edel} result for set-valued uniformly local contraction \cite{nadler}. Later, Sultana and Vetrivel \cite{mam} extended the Theorem \ref{mizo_thm} due to Mizoguchi and Takahashi \cite{mizo} for set-valued uniformly local contraction. The authors \cite{mam} derived the succeeding theorem.
\begin{theorem}\cite{mam}\label{mizo_local_thm}
	Suppose that $(Z,d)$ is complete $r$-chainable and a set-valued map $F:Z\to CB(Z)$ in order that for each $y,z\in Z$ having $d(y,z)<r$ (where $r>0$) fulfils,
	\begin{equation}\label{mizo_local}
		H(Fy,Fz)\leq k(d(y,z))d(y,z),
	\end{equation}
	where $k\in W$. Then there is an element $z^*\in Fz^*$.
\end{theorem}

It is worth to note that every Mizoguchi-Takahashi local contraction \cite{mam} is not a Nadler \cite{nadler} local contraction, see \cite[Example 3.1]{mam_equivalence}. But, Sultana and Qin \cite{mam_equivalence} demonstrated that these two contractions are equivalent, whenever the metric space $(Z,d)$ is metrically convex complete in the year 2019.

On the other hand, Feng and Liu \cite{feng} improved Nadler's \cite{nadler} theorem for the mappings $F$ from $Z$ into $CL(Z)$ in the year of 2006, which is stated below.
\begin{theorem}\cite{feng}\label{feng_thm}
	Suppose that $(Z,d)$ is a complete metric space and a set-valued map $F:Z\to CL(Z)$ in order that for each $y\in Z$, there is $z\in I_{c\in (0,1)}^y=\{x\in Fy:cd(x,y)\leq D(y,Fy)\}$,
	\begin{equation}\label{feng_eq}
		D(z,Fz)\leq \alpha d(y,z)~\textrm{where}~ 0\leq \alpha\leq c<1.
	\end{equation}
	Then there an element $z^*\in Fz^*$ if $h:Z\to \mathbb{R}$ by $h(z)=D(z,Fz)$ is lower semicontinuous.
\end{theorem}

In this manuscript, inspired by Eldred et. al \cite{eldred} and Sultana-Qin \cite{mam_equivalence}, we establish that any set-valued contraction due to Feng-Liu \cite{feng} follows from the Mizoguchi-Takahashi local contraction \cite{mam_equivalence} on a metrically convex metric space. On any arbitrary metric spaces, this result need not be true, an example is provided for this purpose. The fixed points for the set-valued mappings fulfil the equation $(\ref{mizo_local})$ are also studied on metrically convex spaces. Furthermore, the equivalence of Mizoguchi-Takahashi local contraction \cite{mam_equivalence} and Nadler local contraction \cite{nadler} is established in a compact metric space. An invariant approximation result is also discussed through our main theorem.
\section{Notations and Definitions}
This segment contains with some mathematical definitions and notations, which are needful up to the end of this paper. For the metric space $(Z,d)$, let $u,v$ be any two points in $Z$. Then a point $w$ in $Z$, where $u\neq w\neq v$, is called metrically \cite{eldred} between $u$ and $v$ if $d(u,v)=d(u,w)+d(w,v)$. The space $(Z,d)$ is called metrically convex \cite{eldred} if for every pair of elements $u$ and $v$ in $Z$, there is an element $w\in Z$ in order that $w$ is metrically between $u$ and $v$. The below stated lemma is regarding a property of metrically convex space, which is needful for establish our main theorem.

\begin{lemma}\label{metrically}\cite{eldred}
	Assume that $(Z,d)$ is metrically convex complete and $u,v\in Z$ are any arbitrary points. Then there is $[a_1,a_2]\subset \mathbb{R}$ and an isomorphism $\varphi:[a_1,a_2]\to Z$ in order that $\varphi(a_1)=u$, $\varphi(a_2)=v$.
\end{lemma}

Again for the metric space $(Z,d)$, consider a self map $\Psi$ on $Z$. Then $\Psi$ is defined as uniformly local contractions \cite{edel} if for each $y,z \in Z$ having $d(y,z)<r$ fulfils $d(\Psi (y),\Psi (z))\leq \alpha d(y,z)$, where $0\leq \alpha <1$ and $r>0$. Furthermore, for $y,z\in Z$, a $r$-chain \cite{edel} between $y$ and $z$ means, there is a finite sequence $y=z_0,z_1,\cdots,z_n=z$ in $Z$, such that for every $j\in\{1,2,\cdots,n\}$, $d(z_{j-1},z_j)<r$. The space $(Z,d)$ is stated as $r$-chainable \cite{edel} whenever there is a $r$-chain between every points $y$ and $z$ in $Z$.

\section{Main Result}
Now we commence our main theorem about that any Mizoguchi-Takahashi local contraction \cite{mam}, that is, the mappings that meet the equation (\ref{mizo_local}), is actually a Feng-Liu contraction \cite{feng} on a metrically convex metric space $(Z,d)$.
\begin{theorem}\label{main}
	Suppose that $(Z,d)$ is a complete metric space having metrically convexity property and a set-valued map $F:Z\to CB(Z)$ in order that for $y,z\in Z$ having $d(y,z)<r$ (where $r>0$),
	\begin{equation}\label{MT}
		H(Fy,Fz)\leq k(d(y,z))d(y,z)
	\end{equation}
	where $k\in W$. Then $F$ becomes a Feng-Liu contraction, that is, for each $y^{'}\in Z$, there is $z^{'}\in I_c^{y^{'}}$ fulfils,
	\begin{equation*}
		D(z^{'},Fz^{'})\leq \alpha d(y^{'},z^{'})~\textrm{where}~0\leq \alpha <c<1.
	\end{equation*}
\end{theorem}
\begin{proof}
	Our main goal is to show that there is certain $\alpha$, where $0\leq \alpha<c<1$ in order that for $y\in Z$, $z\in I_c^y$ meets $D(z,Fz)\leq \alpha d(y,z)$. Suppose that
	\begin{equation}
		M=\{m>0: \sup\{k(d(y,z)):d(y,z)\in [0,m]\}=1\}.
	\end{equation}
	Then there are two possibilities that could happen. Firstly, we assume that the set $M$ is void. Then for each positive $m$, the supremum value of $k(d(y,z))$ is strictly less than $1$, where $d(y,z)\in [0,m]$.
	Since $r>0$ is fixed, we assume that
	\begin{equation}\label{set}
		\sup\{k(d(y,z)):d(y,z)\in [0,r]\}=\alpha<1.
	\end{equation}
	Evidently, $\alpha\in [0,1)$. Let $y^{'}\in Z$ be a fixed element. Then there is an element $z^{'}\in I_c^y{^{'}}$ where $c>\alpha$, due to the fact that $F(y^{'})\subset CL(Z)$. Consequently, for chosen $y^{'}$ and $z^{'}$, an isometry $\varphi_1:[t_1,t_2]\to Z$ exists satisfies $\varphi_1(t_1)=y^{'}$ and $\varphi_1(t_2)=z^{'}$ (where $[t_1,t_2]\subset \mathbb{R}$ ) according to the Lemma \ref{metrically}.
	Now for a fixed positive real number $q<r$, we can always find $L\in \mathbb{N}$ satisfying $t_1+Lq<t_2$ and $t_2\leq t_1+(L+1)q$. Consequently,
	\begin{equation}\label{isomorphism}
		\begin{aligned}
			d(y^{'},z^{'})&=d(\varphi_1(t_1),\varphi_1(t_2))=|t_1-t_2|\nonumber\\
			&=|t_1-(t_1+q)|+|(t_1+q)-(t_1+2q)|+\cdots+|(t_1+Lq)-t_2|\nonumber\\
			&=d(y^{'},\varphi_1(t_1+q))+d(\varphi_1(t_1+q),\varphi_1(t_1+2q))+\cdots+\nonumber\\
			&\qquad d(\varphi_1(t_1+Lq),z^{'}).
		\end{aligned}
	\end{equation}

	Again, it is simple to observe that $d(z^{'},\varphi_1(t_1+Lq))\leq q<r$ and $d(\varphi_1(t_1+nq),\varphi_1(t_1+(n+1)q))=q<r$ for every $0\leq n<L$, that is every terms of right hand side of the equation (\ref{isomorphism}) is less than $r$. Therefore for chosen $y^{'}\in Z$ and $z^{'}\in I_c^{y^{'}}$, it yields

		\begin{align*}
			D(z^{'},Fz^{'})&\leq H(Fy^{'},Fz^{'})\\
			&\leq H(Fy^{'},F(\varphi_1(t_1+q)))+\cdots+H(F(\varphi_1(t_1+Lq)),Fz^{'})\\
			&\leq  \alpha[d(y^{'},\varphi_1(t_1+q))+\cdots+d(\varphi_1(t_1+Lq),z^{'})]\\

			&\leq \alpha d(y^{'},z^{'}).\qquad[\text{using}~(\ref{isomorphism}) ]
		\end{align*}

	Next, we assume that the set $M$ is non-void and $m_0$ is the infimum of $M$. If possible, let $m_0=0$. Then we are able to find two sequences $\{d(y_n,z_n)\}_n$, $\{k(d(y_n,z_n))\}_n$ in order that $d(y_n,z_n)\to 0$ and $k(d(y_n,z_n))\to 1$, which contradicts to the criteria of the map $k\in W$. Hence $m_0\neq 0$. Then for a chosen positive real $q_0<m_0$, we have
	\begin{equation}\label{set1}
		\sup\{k(d(x_1,x_2)):d(x_1,x_2)\in [0,q_0]\}=\alpha_0<1.
	\end{equation}
	Let $y^{'}\in Z$ be a fixed element. Subsequently, there is $z^{'}\in I_c^{y^{'}}$, where $c>\alpha_0$. Then for that $y^{'}$ and $z^{'}$, there is an isometry $\varphi_2:[s_1,s_2]\to Z$ in order that $\varphi_2(s_1)=y^{'}$ and $\varphi_2(s_2)=z^{'}$, where $[s_1,s_2]\subset \mathbb{R}$. Consequently, for certain fixed positive number $a<\min\{q_0,r\}$, we can find $L^{'}\in \mathbb{N}$ satisfying $s_1+L^{'}a<s_2$ and $s_2\leq s_1+ (L^{'}+1)a$.
	Evidently,
	\begin{equation}\label{isomorphism1}
		\begin{aligned}
			d(y^{'},z^{'})&=d(\varphi_2(s_1),\varphi_2(s_2))=|s_1-s_2|\nonumber\\
			&=|s_1-(s_1+a)|+|(s_1+a)-(s_1+2a)|+\cdots+|(s_1+L^{'}a)-s_2|\nonumber\\
			&=d(y^{'},\varphi_2(s_1+a))+d(\varphi_2(s_1+a),\varphi_2(s_1+2a))+\cdots+\nonumber\\
			&\qquad d(\varphi_2(s_1+L^{'}a),z^{'}).
		\end{aligned}
	\end{equation}
	Again, it is simple to check that $d(z^{'},\varphi_2(s_1+L^{'}a))\leq a$ and $d(\varphi_2(s_1+na),\varphi_2(s_1+(n+1)a))=a$ for every $0\leq n<L^{'}$, that is each terms of right hand side of the equation (\ref{isomorphism1}) is less than or equal to $a$. Therefore for that chosen $y^{'}\in Z$ and $z^{'}\in I_c^{y^{'}}$, it yields

	\begin{align*}
		D(z^{'},Fz^{'})\leq& H(Fy^{'},Fz^{'})\\
		\leq&H(Fy^{'},F(\varphi_2(s_1+a)))+\cdots+H(F(\varphi_2(s_1+L^{'}a)),Fz^{'})\\
		\leq & \alpha_0[d(y^{'},\varphi_2(s_1+a))+\cdots+d(\varphi_2(s_1+L^{'}a),z^{'})]\\

		\leq&\alpha_0 d(y^{'},z^{'})\qquad[\text{using}~(\ref{isomorphism1})].
	\end{align*}
	Thus $F$ becomes a Feng-Liu contraction.
\end{proof}

In the succeeding theorem, we establish the fixed points for Mizoguchi-Takahashi local contraction on a metrically convex space. Indeed, this result follows from our above mentioned Theorem \ref{main} and the Theorem \ref{feng_thm} due to Feng and Liu \cite{feng}. Although, this theorem is established by Sultana and Qin \cite[Corollary 3.1]{mam_equivalence} in the year 2019.
\begin{theorem}\label{cor1}\cite{mam_equivalence}
	Suppose that $(Z,d)$ is metrically convex complete and a set-valued map $F:Z\to CB(Z)$ meets the equation $(\ref{mizo_local})$. Then presence of $z^*\in Z$ in order that $z^*\in Fz^*$ can be assured.
\end{theorem}
\begin{proof}
	Since $F$ meets the equation $(\ref{mizo_local})$ and $Z$ is metrically convex, then by Theorem \ref{main}, the map $F$ fulfils the Feng-Liu contraction. As a consequence, for each $y\in Z$, $z\in I_c^y$ where $0\leq \alpha<c<1$ fulfils
	\begin{eqnarray*}
		D(z,Fz)\leq \alpha d(y,z).
	\end{eqnarray*}
	Now our aim is to establish that $h:Z\to \mathbb{R}$ by $h(z)=D(z,Fz)$ is lower semicontinuous. Let us assume $\{q_n\}_n\in Z$ in order that $q_n\to q$. Then we achieve $l\in \mathbb{N}$ in order that $d(q_n,q)<r$, for each $n\geq l$. Now for every $n\geq 1$ and $z\in Fq_n$,
	\begin{eqnarray}
		D(q,Fq)&\leq&d(q,q_n)+d(q_n,z)+D(z,Fq)\nonumber\\
		&\leq&d(q,q_n)+d(q_n,z)+H(Fq_n,Fq)\nonumber\\
		&\leq&d(q,q_n)+D(q_n,Fq_n)+H(Fq_n,Fq).
	\end{eqnarray}
	As for all $n\geq l$, $d(q_n,q)<r$, then the last inequation leads to
	\begin{eqnarray*}
		D(q,Fq)&\leq& d(q,q_n)+D(q_n,Fq_n)+k(d(q_n,q))d(q_n,q)\\
		&< & d(q,q_n)+D(q_n,Fq_n)+d(q_n,q).
	\end{eqnarray*}
	Taking $n\to \infty$, we conclude that $D(q,Fq)\leq\liminf_{n\to \infty}D(q_n,Fq_n)$. Now applying the Theorem \ref{feng_thm}, we are able to find $z^*\in Z$ in order that $z^*\in Fz^*$.
\end{proof}

In 2009, Eldred et. al \cite[Theorem 2.4]{eldred} also established that Mizoguchi-Takahashi contraction \cite{mizo} follows Nadler's contraction \cite{nadler}, whenever the metric space is compact. In the below stated theorem, we establish the equivalence between the Mizoguchi-Takahashi local contraction \cite{mam} and the Nadler local contraction \cite{nadler} in a metric space which is compact.
\begin{theorem}
	Let $(Z,d)$ be compact and a set-valued map $F:Z\to CB(Z)$ meets the equation $(\ref{mizo_local})$. Then for each $y,z\in Z$ having $d(y,z)<r$ (where $r>0$),
	\begin{equation*}
		H(Fy,Fz)\leq \alpha\, d(y,z)\quad \text{where}~~0
		\leq\alpha< 1.
	\end{equation*}
\end{theorem}
\begin{proof}
	Let $0<q<r$ be a fixed real number. Now consider a real number $m$, which is the supremum of $\frac{H(Fy,Fz)}{d(y,z)}$, where the elements $y$ and $z$ lies in $Z$ having $0< d(y,z)\leq q$.

	We want to show that $m<1$. If possible let $m=1$, then there is $\{(y_n,z_n)\}_n\in Z\times Z$ having $d(y_n,z_n)\in (0,q]$ fulfils
	\begin{equation}\label{limit}
		\lim_{n\to\infty}\frac{H(Fy_n,Fz_n)}{d(y_n,z_n)}=1.
	\end{equation}
	On the account of $d(y_n,z_n)\leq q<r$, we have $\frac{H(Fy_n,Fz_n)}{d(y_n,z_n)}\leq k(d(y_n,z_n))<1$. Taking $n\to \infty$, it yields that $k(d(y_n,z_n))\to 1$. If $d(y_n,z_n)$ converges to $0$, then it contradicts the criteria of $k\in W$. Hence we assert that $\{d(y_n,z_n)\}_n$ converges to a positive real number.

	As $\{(y_n,z_n)\}_n\in Z\times Z$ and $Z$ is compact, then there is $\{(y_{n_l},z_{n_l})\}_{l\in \mathbb{N}}$ in order that $(y_{n_l},z_{n_l})\to (y^*,z^*)$. Consequently, $d(y_{n_l},z_{n_l})\to d(y^*,z^*)$. Since $d(y_n,z_n)\in (0,q]$, then $d(y^*,z^*)\in (0,q]$, that is $d(y^*,z^*)\leq q<r$. Now
	\begin{eqnarray*}
		\lim_{l\to\infty}\frac{H(Fy_{n_l},Fz_{n_l})}{d(y_{n_l},z_{n_l})}&=&\frac{H(Fy^*,Fz^*)}{d(y^*,z^*)}\\
		&\leq & \frac{k(d(y^*,z^*))d(y^*,z^*)}{d(y^*,z^*)}\qquad [\because~d(y^*,z^*)<r]\\
		&=&k(d(y^*,z^*))<1.
	\end{eqnarray*}
	This contradicts the equation $(\ref{limit})$. Thus we can conclude that $m<1$, that is there is an non-negative $\alpha$ in order that $m=\alpha<1$. Therefore for each $y,z\in Z$ having $d(y,z)<r$, we have
	\begin{equation*}
		H(Fy,Fz)\leq k(d(y,z))d(y,z)\leq \alpha d(y,z)\quad\text{where}~\alpha\in [0,1).
	\end{equation*}
\end{proof}

Many mathematicians are interested to deal the invariant best approximation problems through fixed point, see (\cite{thagafi,pathak}). In \cite{kamran1,kumam}, the authors established invariant best approximation theorems for set-valued mappings. In the following, using the Theorem \ref{feng_thm} due to Feng-Liu \cite{feng} and our main Theorem \ref{main}, we derive a invariant best approximation theorem for a set-valued mapping which meets the criteria $(\ref{mizo_local})$ on a metrically convex space. Let $P$ be a non-void subset of a normed linear space $(Z,\|.\|)$ and $p\in Z$. Then a collection of elements $y\in P$ in order that $D(p,P)=\|y-p\|$ is described as $B_P(p)$ and it is stated as a best $P$-approximates of $p$ over $P$.

\begin{theorem}
	Let $P$ be a subset of a metrically convex normed linear space $(Z,\|.\|)$ and a set-valued map $F:P\to CB(P)$ in order that for every $y,z\in P$ having $d(y,z)<r$ fulfils $(\ref{mizo_local})$. Suppose that the below stated criteria hold:
	\begin{enumerate}
		\item [(i)] the set $B_P(p)$ is a complete subset of $P$,
		\item [(ii)] for every $x\in B_P(p)$, $\sup_{y\in Fx}\|y-p\| \leq \|x-p\|$.
	\end{enumerate}
	Then there is $z^*\in B_P(p)\cap \{z\in Z: z\in Fz\}$ if the function $h:B_P(p)\to \mathbb{R}$ is lower semicontinuous.
\end{theorem}
\begin{proof}
	As $Z$ is metrically convex and the set-valued map $F:P\to CB(P)$ fulfils $(\ref{mizo_local})$ for each $y,z\in P$ having $d(y,z)<r$, then by the Theorem \ref{main}, for each $y\in B_P(p)\subseteq P$, $z\in I_c^y$ meets $D(z,Fz)\leq \alpha d(y,z)$, where $0\leq \alpha<c<1$. Now let $x\in B_P(p)$ and $z\in Fx$. As $x\in B_P(p)$, then $D(p,P)=\|p-x\|$. Again from $(ii)$ we have
	\begin{equation*}
		\|z-p\|\leq \sup_{y\in Fx}\|y-p\| \leq \|x-p\|=D(p,P).
	\end{equation*}
	Thus we can conclude that $\|z-p\|=D(p,P)$. Therefore $z\in B_P(p)$. Hence $Fx\subseteq  B_P(p)$ for each $x\in B_P(p)$, that is $F( B_P(p))\subseteq  B_P(p)$. Immediately it observe that $Fx$ are also closed for any $x\in B_P(p)$, on the account of $Fx\subset CL(P)$, for each $x\in P$. Consequently $F|_{ B_P(p)}: B_P(p)\to CL( B_P(p))$. Thus $F$ fulfils all the criteria of the Theorem \ref{feng_thm} on the set $B_P(p)$ and hence
	$$w^*\in \{z\in B_P(p):z\in F|_{ B_P(p)}z\}= \{z\in Z: z\in Fz\}\cap B_P(p).$$
\end{proof}

The succeeding example illustrates that any uniformly local contraction (\ref{mizo_local}) due to Sultana and Vetrivel \cite{mam} does not follow  Feng-Liu contraction \cite{feng} on any arbitrary metric spaces.

\begin{example}
	Consider $Z=[0,\frac{1}{2}]\cup \{n\in \mathbb{N}: n \geq 2\}$ and $d:Z\times Z\to \mathbb{R}$ is defined by
	\begin{center}
		$d(y,z)=
		\begin{cases}
			~0 &\mbox{if}~y=z,\\
			~|y-z| &\mbox{if}~0\leq y,z\leq \frac{1}{2},\\
			~y+z &\mbox{if one of}~ y,z\notin[0,\frac{1}{2}].
		\end{cases}$
	\end{center}
	Indeed $(Z,d)$ is not a metrically convex metric space. Now define a map $F:Z\to CB(Z)\subseteq CL(Z) $ in order that
	\begin{center}
		$Fz=
		\begin{cases}
			~\left\{\frac{1}{2}z^2 \right\} &\mbox{if}~z\in[0,\frac{1}{2}],\\
			~\{7,2z-1\} &\mbox{if}~ z\in\{2,3,4,\cdots\}.
		\end{cases}$
	\end{center}
	Our goal is to show that $F$ is not a Feng-Liu contraction but a Mizoguchi-Takahashi local contraction. Now for $3\in Z$, $F3=\{7,5\}$. Then there are two cases that arise.
	\begin{enumerate}
		\item [(i)] Choose $7\in \{7,5\}$, then $F7=\{7,13\}$. Now we observe that $D(7,F7)=\inf\{d(7,7),d(7,13)\}=\inf\{14,20\}=14$. Again $d(3,7)=10$. Therefore $$D(7,F7)> d(3,7).$$

		\item[(ii)] Choose $5\in \{7,5\}$, then $F5=\{7,9\}$. Now we see that $D(5,F5)=\inf\{d(5,7),d(5,9)\}=\inf\{12,14\}=12$. Again $d(3,5)=8$. Therefore $$D(5,F5)> d(3,5).$$
	\end{enumerate}
	Hence for $y=3\in Z$ and every $z\in F3$, we obtain that
	$$D(z,Fz)>d(y,z).$$
	This indicates that $F$ is not a Feng-Liu contraction. On the other side, when $d(y,z)<\frac{1}{2}$, then it is easy to visualize that $y,z\in [0,\frac{1}{2}]$. Now
	\begin{eqnarray*}
		H(Fy,Fz)&=&H\left(\left\{\frac{1}{2}y^2\right\},\left\{\frac{1}{2}z^2\right\}\right)\\
		&\leq& \frac{1}{2}(y+z)d(y,z).
	\end{eqnarray*}
	Consider a map $k:[0,\infty)\to [0,1)$ by
	\begin{center}
		$k(s)=
		\begin{cases}
			~\frac{3}{4 } &\mbox{if}~s\in[0,1/2)\\
			~0 &\mbox{if}~ s\in[1/2,\infty).
		\end{cases}$
	\end{center}
	Therefore $k\in W$ and hence for $d(y,z)<\frac{1}{2}$,
	$$H(Fy,Fz)\leq k(d(y,z))d(y,z).$$
\end{example}

\section*{Acknowledgment}

The first author would like to acknowledge the Ministry of Human Resource Development, India for providing financial assistance during the research work.


\section{Introduction}

A crucial process in printed circuit board assembly is the visual inspection of electronic components for potential defects. This can help avoid functional failure of devices, user data leakage, or even system control taken by adversaries ~\citet{cryptoeprint:2020/366}. Given that there can be hundreds of electronic components on a given PCB, the task of visual inspection can be extremely time-consuming and prone to operator error, especially during large assembly runs. Therefore, the ability to automatically detect different electronic components on a PCB board for automated inspection purposes is highly desired. As a result, there has been significant interest in the research community in automatic PCB component detection, particularly leveraging deep learning ~\citet{kuo2019data, electronics11081183}. However, one consideration that has been largely left unexplored in research literature in this area is computational efficiency, which is particularly critical for real-world visual quality inspection scenarios involving high-volume, high-throughput electronics manufacturing use-cases under constrained edge computing resources.

Motivated by the need for both high efficiency and high accuracy for automatic PCB component detection, this study explores efficient deep neural network object detection architectures for the purpose of automatic PCB component detection on the edge.  As a result of this exploration, we introduce PCBDet, a highly efficient, performant self-attention deep neural network architecture design. This architecture notably makes use of the recently introduced AttendNeXt backbone, which has been shown to achieve state-of-the-art performance for TinyML, and is integrated here into RetinaNet~\citet{attendnext, retinanet}.

The paper is organized as follows. In Section 2, details about the architecture of PCBDet, the training procedure, the evaluation procedure, and the data used for training and evaluation are described. In Section 3, experimental results in terms of component detection performance, model size, and inference speed on different computing hardware are presented. Finally, conclusions are drawn and future directions are discussed in Section 4.

\section{Methods}

\begin{figure}[t]
        \makebox[\textwidth][c]
		{\includegraphics[scale = 0.37]{PCBDet.jpg}}%
	\caption{Overview of PCBDet architecture. PCBDet consists of A) a double-condensing attention condenser feature encoder feeding a FPN, and B) classification and box regression convolutional sub-nets for bounding box prediction. }
	\label{fig1}
\end{figure}

\subsection{Dataset}

To appropriately explore the quality and impact of our network design, a dataset that can facilitate the training and validation of robust models is essential. With this in mind, models were trained and tested using the DSLR images of the FICS-PCB dataset, a public, comprehensive, and diverse PCB component dataset that includes a number of challenging cases. The dataset itself consists of a total of 31 PCB boards containing over 77 thousand components to detect, with capacitors and resistors being the most widely represented pieces~\citet{cryptoeprint:2020/366}. Each board in the FICS-PCB dataset was further truncated into square patches, reducing train-time resource demand but maintaining component-wise resolution, with each patch being reshaped for additional size reduction further in the input pipeline. Figure 2 demonstrates several examples of such patches with annotated ground truth component labels.

Division of the dataset into distinct train, validation, and test splits is another crucial element in confirming the soundness of our experimentation. The preservation of exclusivity in the test set here is integral, since it allows for performance evaluation on a strict holdout set, and thus all of the image patches extracted from seven of the 31 PCBs were used as the test set. With the exception of one of the boards, which was excluded due to a lack of DSLR pictures, patches from the remaining boards (which total to 23) were used for the train/validation splits. From this set, 87.5\% of the patches were taken for the train set, while the other 12.5\% were used for the validation set (for post-epoch performance validation).

\begin{figure}[t]
	\makebox[\textwidth][c]
		{\includegraphics[scale = 0.32]{6_pcbs.jpg}}%
	\caption{Examples of PCB patches annotated with ground truth bounding boxes.}
	\label{fig2}
\end{figure}

\subsection{Architecture}

As seen in Figure 1, the proposed PCBDet possesses an efficient self-attention architecture design inspired by two different architecture design paradigms: 1) RetinaNet bounding box prediction structure~\citet{retinanet}, and 2) double-condensing attention condenser architecture~\citet{attendnext}. As a single-stage detector architecture design, the RetinaNet structure encompasses a more efficient object detection process when compared to state-of-the-art two-stage object detectors like R-CNN~\citet{retinanet}. RetinaNet has also seen increased performance when compared to one-stage detectors such as SSD or YOLO~\citet{retinanet}. As such, the proposed PCBDet architectural design takes inspiration from the RetinaNet structure, looking to adopt an efficient single-stage approach without seeing substantial tradeoffs in performance.

Without an efficient backbone, however, our network cannot maximize on the possible efficiency-based benefits of the RetinaNet framework. More specifically, the backbone architecture within a RetinaNet structure is the feature encoder that feeds into the convolutional sub-nets, and while a larger, complex backbone may enable increased performance gains, this can lead to substantial losses in efficiency. The design of a small, efficient backbone architecture is thus crucial in creating an effective and efficient object detection network.

PCBDet's backbone architecture takes inspiration from the AttendNeXt double-condensing attention condenser architecture design, which has shown top of the line performance among other state-of-the-art efficient architectures for ImageNet classification~\citet{attendnext}. This self-attention architecture design features double-condensing attention condenser modules, applied within a convolutional network structure, to increase the speed and efficiency of standard convolutional architectures for feature extraction, thus serving as the basis for an efficient backbone for RetinaNet~\citet{attendnext}. The AttendNeXt feature encoder used in our study was first pretrained on ImageNet, establishing a basis for the weights to be used in our object detection task. The classification head was then removed, and stage outputs were taken as inputs for a feature pyramid network (FPN), whose layers respectively feed into classification and regression subnets as per general RetinaNet structure. To further increase the efficiency of our designed network, the first of four stages of the feature encoder was omitted from the construction of the FPN, decreasing the amount of subnet operations performed per pass. The resultant network from this design process is dubbed PCBDet.

We also tried integrating an EfficientNetB4-based backbone into RetinaNet and compared the performance of this architecture with PCBDet. EfficientNets are a family of convolutional models generally designed for, as the name implies, efficiency, with the model seeing upscaling as it progresses from B0 through B7~\citet{efficientnet}. EfficientNet-B4 shows improved top-1 ImageNet performance when compared to other state-of-the-art convolutional classifiers while maintaining a lower number of parameters, and was thus chosen as an efficient but potent backbone to explore with RetinaNet~\citet{efficientnet}. As with PCBDet's AttendNeXt backbone, the EfficientNet feature encoder had its classification head removed and block outputs were used as inputs for an FPN. To provide a fair point of comparison for PCBDet, the integration of this feature encoder with the FPN once again seeks to achieve greater efficiency, with the first three of seven blocks of the EfficientNet-B4 feature encoder being excluded from FPN construction. The network designed here is referred to as EfficientNet-Det.

\subsection{Training}

Proper exploration of our architectures requires thorough training, and for compact networks such as PCBDet in particular, slower, gradual weight learning is crucial to appropriately search for effective weights in the training process.
As such, PCBDet and EfficientNet-Det were each trained for 300 epochs, with a base learning rate of 2e-4 and a proprietary learning rate scheduler, along with Adam optimization. While potential overfitting could arise from slow, gradual learning, this issue was combated with the use of image augmentation, including vertical and horizontal flipping and translation, colour degeneration, and random cutouts (patch removal), as well as the monitoring of network performance on the validation set. It is also essential to address the disproportionate representation of components in the FICS-PCB dataset. To do so, network training uses the focal loss metric, which accounts for class imbalances by adding a focusing parameter to the standard cross-entropy loss, resulting in greatly decreased weighting for easy, well-classified data points~\citet{retinanet}.

During training, the first and second blocks of the AttendNeXt feature encoder in PCBDet were frozen, allowing for the encoder to retain its memories of low-level features from ImageNet pretraining while also tuning higher-level feature blocks to better recognize the shapes and objects associated with PCB components. Similarly, the first four of seven blocks were frozen for the EfficientNet-B4 encoder in EfficientNet-Det.

\subsection{Evaluation}
Given the goal of efficient model design, we need a method that can effectively measure the complexity and compactness of a model. Inference time, or the time taken per forward pass, is a method that can reveal how quickly a network can perform as a predictor; in our work, inference time was taken for both PCBDet and EfficientNet-Det using an Intel Core i7-10750H processor and a NVIDIA GeForce RTX 1650 Ti, both within a Dell XPS 15 9500 laptop, as well as a Jetson Nano and a 64-bit ARM Cortex A72 processor, altogether providing an image of the on-the-edge inference speeds of the two networks. The number of parameters was also taken for each of the two networks, providing an additional measure for model compactness.

The predictive power for bounding boxes of our networks is another necessary measure to analyze, in order to compare the model performance of the PCBDet and EfficientNet-Det architectures. This performance assessment was realized using the mean average precision (mAP) for bounding box predictions, over IOU thresholds from 0.5 to 0.95 with a 0.05 step size, commonly known as mAP@[0.5:0.95]; this mAP metric is also known as COCO mAP, the standard performance metric for COCO challenges~\citet{coco}.
Averaging performance over a range of IOU thresholds provides a more generalized sense of object detection performance across resolutions, as lower IOU thresholds test for more roughly correct box predictions while higher thresholds solely reward exact bounding box location. The validation and test performances of PCBDet and EfficientNet-Det were taken to be the mAP@[0.5:0.95] on the validation and holdout test sets respectively; this general AP metric helps to determine the predictive accuracy of our networks during and after training.

While individual differences in network performance and compactness can be seen through inference time and COCO mAP measures, a collective assessment can provide a better picture of the difference in the accuracy-complexity balance achieved in PCBDet and EfficientNet-Det. This unified analysis can be performed using the NetScore metric, which acts as a quantitative method of assessing this very balance~\citet{netscore}. In this calculation, the coefficient values used were $\alpha$ = 2, $\beta$ = 1, and $\gamma$ = 1, in accordance with the original NetScore study~\citet{netscore}. The inference-time multiply-accumulate (MAC) operations measure was also replaced with the experimental inference time (in seconds) using the ARM Cortex A72 in the calculation, as this experimental metric of complexity provides a more practical measure for edge performance than MAC operations, while the COCO mAP was used as the accuracy metric for the calculation. The calculation used for NetScore was

\[NetScore = \frac{(mAP*100)^{2}}{(MParams)(Inference \; \, time \;  (s))}\]

\section{Results}
The efficacy of the proposed PCBDet for PCB component detection is compared here with EfficientNet-Det across the following metrics: 1) COCO mAP, 2) model size, and 3) inference speed on various low-power computing hardware.

\begin{figure}[H]
	\begin{center}
		\includegraphics[scale = 0.8]{pcbdet_map.jpg}

	\end{center}
	\caption{COCO mAP on validation and test data for PCBDet and EfficientNet-Det.}
	\label{fig3}
\end{figure}

\noindent \textbf{COCO mAP}. It can be observed in Figure 3 that the proposed PCBDet achieves noticeable gains in terms of test and validation COCO mAP by approximately 4 \% and 2 \%, respectively,  when compared to EfficientNet-Det.   This gain in mAP is particularly interesting especially given the fact that PCBDet is significantly smaller and faster than EfficientNet-Det, which we will discuss in greater detail. As such, these results illustrate that a high level of performance can be achieved with the proposed PCBDet for the purpose of automatic PCB component detection on the edge.

\noindent \textbf{Model Size}.  As shown in Figure 4, it can be observed that the proposed PCBDet possesses less than half the number of total parameters, as well as trainable parameters, when compared to EfficientNet-Det. This is particularly important for edge based applications such as automatic on-the-edge PCB component detection where memory resources are limited.

\noindent \textbf{Inference Speed}. As shown in Figure 5, it can be observed that the proposed PCBDet is more than 30\% faster than EfficientNet-Det on the NVIDIA Geforce RTX 1650 Ti, with an even greater speed gain on slower hardware such as the Jetson Nano (almost 65\% faster) and Intel Core i7-10750H (over 45\% faster).  As seen in  Figure 6, the performance gains of the proposed PCBDet on lower-power hardware were especially apparent when evaluated on an ARM Cortex A72, where PCBDet was more than 2$\times$ faster than EfficientNet-Det. These inference speed results demonstrate the efficacy of the proposed PCBDet for high-throughput PCB component detection on the edge.

Finally, using the above results, PCBDet was found to achieve a NetScore of 28.2670 while EfficientNet-Det achieved a NetScore of 13.5749, supporting our findings that PCBDet achieves a superior accuracy-complexity balance when compared to EfficientNet-Det.

These results demonstrate overall that PCBDet shows significantly greater efficacy than RetinaNet with EfficientNet-B4, which is currently considered a state-of-the-art backbone for TinyML. Ultimately, we have developed a model for PCB object detection that shows very strong performance despite its small size and high inference throughput.

\begin{figure}[H]
	\begin{center}
		\includegraphics[scale = 0.55]{params.PNG}

	\end{center}
	\caption{Number of total and trainable parameters for PCBDet and EfficientNet-Det.}
	\label{fig4}
\end{figure}

\bigskip

\begin{figure}[H]
	\begin{center}
		\includegraphics[scale = 0.5]{inf_time_1.PNG}

	\end{center}
	\caption{Inference time for PCBDet and EfficientNet-Det across different hardware.}
	\label{fig5}
\end{figure}

\section{Conclusion}
Here, we conducted an exploration of efficient deep neural network object detection architectures for the purpose of automatic PCB component detection on the edge. The resulting network architecture, which we coin PCBDet, methodically integrates the recently introduced AttendNeXt backbone into RetinaNet. This  results in an architecture which can achieve up to a 2x inference speed-up on low power hardware compared to other state-of-the-art efficient architectures, while still achieving a higher mAP on the FICS-PCB benchmark dataset. This makes PCBDet very well-suited for component detection in high-throughput manufacturing scenarios with limited computational resources. Future work may include seeing if a similar strategy involving the methodical use of the AttendNeXt backbone could be employed to develop high performance, efficient deep neural network object detection architectures for other applications.

\noindent
\begin{figure}[t]
	\begin{center}
		\includegraphics[scale = 0.4]{inf_time_2.PNG}

	\end{center}
	\caption{Inference time for PCBDet and EfficientNet-Det on ARM Cortex A72 processor.}
	\label{fig6}
\end{figure}

\bibliographystyle{plainnat}


\section{Introduction}

	%

	%

	%

	%

In recent years, commercial enterprises such as SpaceX, Amazon and Oneweb are building up mega-constellations with hundreds or thousands of low earth orbit (LEO)  satellites to provide global coverage at low-latency and high-bandwidth. As reported by ITU \cite{ITU2021}, nearly $2.9$ billion people (around $37\%$ of the world's population) still do not have access to Internet. More than $70\%$ of the surface of the earth has no terrestrial network coverage (e.g., in ocean, dessert and rain forest areas). Satellite networks (SNs) can not only supplement current terrestrial networks in areas short of Internet infrastructure, but also help provide data communication services for areas where the terrestrial networks are deployed but highly loaded.
Enhanced by inter-satellite links (ISLs) and on-board computing resources, SNs have potential to offer real time applications for ground users with a latency in tens of milliseconds \cite{lai2022spacertc}.

Traditional satellites are customized for particular tasks, without cooperation and resulting in resource under utilization and high operation cost \cite{sheng2017toward}.
It is not flexible for customized SNs to update software or develop new functions, which brings great challenges to rapid adoption of advanced strategies to improve network performance. To address these challenges, virtualization technologies, including software defined networking (SDN) and network function virtualization (NFV), have been introduced into the SNs. In particular, NFV decouples network functions from the dedicated hardwares, thereby, network functions can be virtualized into software components or abstracted as containers, which are referred to as virtual functions (VFs).
In this way, VFs can be flexibly deployed on different satellites to establish customized virtual networks for diversified requested services. Each requested service can specify its targeted VF and quality of service (QoS) requirements such as transmission capacity and end to end delay.
Both the academia and industry efforts have been made on the feasibility of using NFV in the SNs.
In \cite{xu2018software}, NFV is exploited to facilitate the incorporation of new applications in the SNs. The benefits of applying NFV into the SNs are verified through various use cases in \cite{bertaux2015software}.
The implementation of NFV in SNs is investigated in \cite{zhou2019bidirectional}.
The key technologies of NFV have already been verified by a couple of satellites in space, such as Tianzhi 1 and Eutelsat Quantum.
However, research on routing for real time applications over SN with VFs is absent, which is the main focus of our work.

Despite the flexibility and benefits introduced by NFV, enabling real time communications in SNs with NFV still faces many challenges. The multimedia-based real time applications, such as video-conferencing, Internet telephony and interactive VR/AR applications, are delay sensitive in nature. To provision a requested real time service, the routing strategy must have fast response time and satisfy both the QoS and network function requirements under the time-varying topologies.
Therefore, the routing strategy becomes more complex than traditional ones designed for terrestrial networks.
Additionally, there are few efforts to investigate routing in SNs with NFV. In \cite{wang2020sfc, jia2021vnf, yang2021maximum}, the optimization of NFV deployment and routing strategy in the SNs was investigated. However, their \textit{store-wait-forward} data transmission mechanism is designed for delay tolerant tasks, which can incur unacceptable latency in the order of minutes or even hours (especially in large scale SNs) and is not applicable for delay sensitive services.

To fill this gap, we investigate the routing strategy for real time applications in the SNs with NFVs. We identify that such a routing problem over the time-varying graph can be formulated as an integer linear programming problem. The branch and bound method for solving it incurs $\mathcal{O}(|\mathcal{L}^{\tau}| \cdot (3 |\mathcal{V}^{\tau}| + |\mathcal{L}^{\tau}|)^{|\mathcal{L}^{\tau}|})$ time complexity, where $|\mathcal{V}|^{\tau}$ is the number of nodes, and $|\mathcal{L}^{\tau}|$ is the number of links in the snapshot graph within time interval $\tau$. By adopting a k-shortest path-based algorithm, the worst-case time complexity becomes $O(|\mathcal{V}^{\tau}|! \cdot |\mathcal{V}^{\tau}|^3)$, but is fast enough in practical cases. However, its solution can be suboptimal, which may result in compromised acceptance ratio. To overcome this, we design an alternative algorithm by exploiting the special structures of the solution space, which can obtain the optimal solution in polynomial time with a low computation complexity of $\mathcal{O}(3|\mathcal{L}^{\tau}| + (2\log{|\mathcal{V}^{\tau}|}+1) |\mathcal{V}^{\tau}|)$ and perform more stable in large scale networks. Simulations conducted on starlink constellation with thousands of satellites corroborate the effectiveness of the proposed algorithm.

	\section{System Model and Problem Formulation}
			\subsection{Satellite network scenario}

	We consider one typical SN composed of satellites and ground terminals (GTs), which are denoted by $\mathcal{S} = \{\mathbb{S}_1, ..., \mathbb{S}_P\}$ and $\mathcal{O} = \{\mathbb{O}_1, ..., \mathbb{O}_K\}$, respectively. Here $P$ and $K$ are the numbers of elements in $\mathcal{S}$ and $\mathcal{O}$, respectively. Each satellite is equipped with computing and communication hardwares. The computing resources are defined into different types of VFs (such as object recognition, data mining, information encryption and signal processing) by taking advantage of SDN, NFV and micro-services techniques, which guarantees that different types of services can be provisioned in the same SN with flexibility. ISLs are supported in the SN, which are intermittent and predictable with satellites' movement.

	Denote the set of functions supported in the SN as $\mathcal{F} = \{f_1, f_2, ..., f_N\}$, where $N = |\mathcal{F}|$. We assume different numbers and types of functions have been deployed on different satellites in the SN. Specifically, each satellite can provide one or multiple functions, which is different from our previous work \cite{yang2021maximum} assuming only one function is supported per satellite. To prevent the functions of a single satellite from being called by excessive applications, the maximum allowable call number of a function $f_l \in \mathcal{F}$ in a satellite $\mathbb{S}_i \in \mathcal{V}$ is set as $w^{f_l}_{\mathbb{S}_i} $, where $w^{f_l}_{\mathbb{S}_i} \in \{0, 1, 2, 3,...\}$. Specifically, if $w^{f_l}_{\mathbb{S}_i} > 0$, satellite $\mathbb{S}_i$ can handle applications requiring $f_l$ function; if $w^{f_l}_{\mathbb{S}_i} = 0$, no more applications with function $f_l$ can be provisioned.

	Given a time horizon $\mathcal{T}$, we use $\mathcal{A} = \{\mathbb{O}_s, \mathbb{O}_d, f_a, C_a, D_a \}$ to denote the real time communication application requiring function $f_a$, where $\mathbb{O}_s$ denotes the source GT, $\mathbb{O}_d$ denotes the destination GT, $C_a$ is the required transmission bandwidth, and $D_a$ is the maximum acceptable end to end delay from $\mathbb{O}_s$ to $\mathbb{O}_d$. In this work, we assume each application requires one function, and the corresponding data processing is conducted on one single satellite with no split. In other words, an end to end feasible path $p$ for $\mathcal{A}$ must contain at least one satellite providing function $f_a$ for data flow processing, while the other satellites within the path $p$ only relay data without providing functions.
    Besides, the transmission bandwidth and end to end delay requirements are also specified by the user applications, however, how long the application will last is unknown.

	\begin{figure}[t!]
		\centering
				\includegraphics[width=70mm]{System_model_for_satellites_slicing_ICC.pdf}\
		\caption{An example of real time communications over satellite networks. }
	\end{figure}

	\textbf{An example.}
	A scenario of an SN supporting three kinds of VFs is shown in \textcolor{black}{Figure 1}, where three VFs are marked by red square, blue triangle, green pentagon, respectively. There are three real time applications, and each of them has a pair of source and sink GT (e.g., \{s1, d1\}), Qos requirements and the specified function. The allocated paths are presented, and the selected satellites providing functions are highlighted in yellow. For \textit{App\_1} and \textit{App\_3}, simple paths bypassing a satellite deployed with the requested VF are allocated. For \textit{App\_2}, due to the scarcity of its specified VF (i.e., the blue triangle), a non-simple path with repeated nodes is provisioned.

    \subsection{Time-varying graph model}

In our system model, we divide the time horizon $\mathcal{T} = [0, T]$ into variable length time windows using the time division method in \cite{wang2022enhancing}. Each time window is denoted by $\tau = [t_s, t_e]$, where $t_s$ and $t_e$ are the start and end time of $\tau$, and $|\tau| = t_e - t_s $ is its length, $0 \leq t_s < t_e \leq T$. We split $\mathcal{T}$ in such a way that in each time window $\tau$, the link connectivities and their attributes (i.e., link delay and transmission capacity) remain unchanged.

We use the snapshot graphs to describe the dynamic evolution of the topology in discrete time windows.
	The snapshot in time interval $\tau = [t_{i-1}, t_i]$ is depicted as a \textcolor{black}{directed} graph $\mathcal{G}^{\tau} = (\mathcal{V}^{\tau}, \mathcal{L}^{\tau})$, where $\mathcal{V}^{\tau} = \mathcal{S} \cup \mathcal{O}$ is the node set including both the satellites and GTs, and  $\mathcal{L}^{\tau} = \{ (\mathbb{V}_i, \mathbb{V}_j) | \text{ node }\mathbb{V}_i \text{ is in the coverage of node } \mathbb{V}_j\}$ is the set of communication links (i.e., opportunities) including both the ISLs and ground satellite links (GSLs) during time interval $\tau$.

	For each satellite $\mathbb{S}_i \in \mathcal{S}$, the maximum allowable numbers of calls of each function are collected into a set $\mathcal{W}_{\mathbb{S}_i} = \{w^{f_1}_{\mathbb{S}_i}, w^{f_2}_{\mathbb{S}_i}, ..., w^{f_N}_{\mathbb{S}_i}\}$, where $|\mathcal{W}_{\mathbb{S}_i}| = N$.
	Each link $(\mathbb{V}_i, \mathbb{V}_j)$ is characterized by its transmission rate $r^{\tau}_{(\mathbb{V}_i, \mathbb{V}_j)}$ (unit: Mbps) and its delay $D^{\tau}_{(\mathbb{V}_i, \mathbb{V}_j)}$ (unit: ms), where $D^{\tau}_{(\mathbb{V}_i, \mathbb{V}_j)}$ is the sum of the propagation delay, transmission delay and queuing delay.

	\section{Problem Formulation}

	\subsection{Basic constrains for path selection}

		We define the binary variable ${x}^{\tau}_{\mathbb{V}_i, \mathbb{V}_j}$ for each link $(\mathbb{V}_i, \mathbb{V}_j) \in \mathcal{L}^{\tau}$, where ${{x}^{\tau}_{\mathbb{V}_i, \mathbb{V}_j} = 1}$ means the allocated path for the application $\mathcal{A}$ at time interval $\tau$ will pass through link $(\mathbb{V}_i, \mathbb{V}_j)$.
		Inherently, there are the following constraints.
	\subsubsection{\text{Source node constraint}} Obviously, for the source GT, the path must come out from one of its adjacent links, namely,
	\begin{equation}
	\sum_{\mathbb{V}_k \in \mathcal{V}^{\tau} - \{\mathbb{V}_s\}} {x}^{\tau}_{\mathbb{O}_s, \mathbb{V}_k} = 1. \label{c_source_node_constraints}
	\end{equation}
	\begin{equation}
	\sum_{\mathbb{V}_k \in \mathcal{V}^{\tau} - \{\mathbb{V}_s\}} {x}^{\tau}_{\mathbb{V}_k, \mathbb{O}_s} = 0. \label{c_source_node_remove_loop_constraints}
	\end{equation}
	\subsubsection{\text{Sink node constraint}} Similarly, for the destination GT, the path must enter one of its adjacent links, namely,
	\begin{equation}
	\sum_{\mathbb{V}_k \in \mathcal{V}^{\tau} - \{\mathbb{O}_d\}} {x}^{\tau}_{\mathbb{V}_k, \mathbb{O}_d} = 1. \label{c_sink_node_constraints}
	\end{equation}
	\begin{equation}
	\sum_{\mathbb{V}_k \in \mathcal{V}^{\tau} - \{\mathbb{O}_d\}} {x}^{\tau}_{\mathbb{O}_d, \mathbb{V}_k} = 0. \label{c_sink_node_remove_loop_constraints}
	\end{equation}
	Note that the purpose of constraint (\ref{c_source_node_remove_loop_constraints}) and (\ref{c_sink_node_remove_loop_constraints}) is to prevent the loops of a path. However, sub-tours may still occur.

	\subsubsection{\text{Relay node constraint}}
	For any relay node, a path can only pass through it no more than twice. That is, if the data is forwarded by it to reach the function node, the processed data will be transferred to the destination node through it again, or through other relay nodes. Thus, at most two of its incoming links can be included within the planned path, namely,
	\begin{equation}
	\sum_{\mathbb{V}_k:(\mathbb{V}_k, \mathbb{V}_{\xi}) \in \mathcal{L}} {x}^{\tau}_{\mathbb{V}_k, \mathbb{V}_{\xi}} \leq 2, \forall \mathbb{V}_{\xi} \in \mathcal{V}^{\tau} -\{\mathbb{O}_s, \mathbb{O}_d\}. \label{c_incoming_node_constraints}
	\end{equation}
	And similarly, for its outgoing links, no more than two among them can be selected, namely,
	\begin{equation}
	\sum_{\mathbb{V}_k:(\mathbb{V}_{\xi},\mathbb{V}_k) \in \mathcal{L}} {x}^{\tau}_{\mathbb{V}_{\xi}, \mathbb{V}_k} \leq 2, \forall \mathbb{V}_{\xi} \in \mathcal{V}^{\tau} -\{\mathbb{O}_s, \mathbb{O}_d\}. \label{c_outgoing_node_constraints}
	\end{equation}

	\subsubsection{\text{Path integrity constraint}} If a path enters one incoming link of a forwarding node, it must also come out from one of its outgoing links. Thus, $ \forall \mathbb{V}_{\xi} \in \mathcal{V}^{\tau} -\{\mathbb{O}_s, \mathbb{O}_d\}$,
	\begin{equation}
	\sum_{\mathbb{V}_k:(\mathbb{V}_k, \mathbb{V}_{\xi}) \in \mathcal{L}^{\tau}} {x}^{\tau}_{\mathbb{O}_k, \mathbb{V}_{\xi}}  = \sum_{\mathbb{V}_k:(\mathbb{V}_{\xi},\mathbb{V}_k) \in \mathcal{L}^{\tau}} {x}^{\tau}_{\mathbb{V}_{\xi}, \mathbb{V}_k}. \label{c_path_integrity_constraints}
	\end{equation}

	\subsubsection{\text{Path capacity constraint}} All of the transmission capacities of the links in the selected path must be larger than the application $\mathcal{A}$'s required capacity $C_a$, namely,
	\begin{equation}
	 M \cdot (1-{x}^{\tau}_{\mathbb{V}_i, \mathbb{V}_j} ) + {x}^{\tau}_{\mathbb{V}_i, \mathbb{V}_j} \cdot {C}^{\tau}_{\mathbb{V}_i, \mathbb{V}_j} \geq C_a, \forall (\mathbb{V}_i, \mathbb{V}_j) \in \mathcal{L}^{\tau},  \label{c_path_capacity_constraints}
	\end{equation}
	where $M$ is a big constant commonly used in integer programming and can be set as the maximum communication capacity in SN. If ${x}^{\tau}_{\mathbb{V}_i, \mathbb{V}_j} = 1$, then $ M \cdot (1-{x}^{\tau}_{\mathbb{V}_i, \mathbb{V}_j}) = 0$, indicating the capacity of the selected link $ (\mathbb{V}_i, \mathbb{V}_j) $ must be greater than $C_a$. On the contrary, if ${x}^{\tau}_{\mathbb{V}_i, \mathbb{V}_j} = 0$, constraint (\ref{c_path_capacity_constraints}) becomes $M \geq C_a$ which automatically holds thus cancels the restriction on link $ (\mathbb{V}_i, \mathbb{V}_j) $'s capacity.
	\subsubsection{\text{Path delay constraint}} The end-to-end delay of the selected path cannot exceed the required delay bound $D_a$. Intuitively, this constraint can be written as follows.
	\begin{equation}
	\sum_{(\mathbb{V}_i, \mathbb{V}_j) \in \mathcal{L}^{\tau}} {x}^{\tau}_{\mathbb{V}_i, \mathbb{V}_j} \cdot {D}^{\tau}_{\mathbb{V}_i, \mathbb{V}_j} \leq D_a. \label{c_path_delay_constraints}
	\end{equation}
	\subsubsection{\text{Virtual function constraint}} The selected path must contain at least one node providing the required function  $f_a$. Formally, such constraint can be written as follows.
	\begin{equation}
    \sum_{\mathbb{V}_\xi \in \mathcal{S}} [(\sum_{\mathbb{V}_k:(\mathbb{V}_k, \mathbb{V}_{\xi}) \in \mathcal{L}^{\tau}} {x}^{\tau}_{\mathbb{V}_k, \mathbb{V}_{\xi}} )\cdot w^{f_a}_{\mathbb{V}_{\xi}} ] \geq 1. \label{c_path_function_constraints}
	\end{equation}

	\subsection{Constraints for eliminating sub-tours}
	To eliminate possible \textit{sub-tours} in the path, we introduce the integer variable ${y}^{\tau}_{\mathbb{V}_i, \mathbb{V}_j}$ for each link $({\mathbb{V}_i, \mathbb{V}_j})$ to represent its link-order in the path, which can be regarded as the \textit{hop number} of the link $({\mathbb{V}_i, \mathbb{V}_j})$ at time slot $\tau$. If the link ${\mathbb{V}_i, \mathbb{V}_j}$ is included in the path, ${y}^{\tau}_{\mathbb{V}_i, \mathbb{V}_j} = h > 0$ indicates that the link $({\mathbb{V}_i, \mathbb{V}_j})$ is the $h-$hop of the path; otherwise, ${y}^{\tau}_{\mathbb{V}_i, \mathbb{V}_j} = 0$.

	\subsubsection{\text{Binding variables}}
	Firstly, if ${x}^{\tau}_{\mathbb{V}_i, \mathbb{V}_j} = 0$, then ${y}^{\tau}_{\mathbb{V}_i, \mathbb{V}_j} = 0$. Otherwise, $1 \leq {y}^{\tau}_{\mathbb{V}_i, \mathbb{V}_j} \leq \sum_{(\mathbb{V}_i, \mathbb{V}_j) \in \mathcal{L}_{\tau}}{x}^{\tau}_{\mathbb{V}_i, \mathbb{V}_j}$, where $\sum_{(\mathbb{V}_i, \mathbb{V}_j) \in \mathcal{L}_{\tau}}{x}^{\tau}_{\mathbb{V}_i, \mathbb{V}_j}$ is the maximum hop number of the planned path. Such a relationship can be expressed as follows.
	\begin{equation}
	 {x}^{\tau}_{\mathbb{V}_i, \mathbb{V}_j} \leq {y}^{\tau}_{\mathbb{V}_i, \mathbb{V}_j} \leq {x}^{\tau}_{\mathbb{V}_i, \mathbb{V}_j} \cdot M. \label{c_bind_xy_constraints}
	\end{equation}
	\begin{equation}
	{y}^{\tau}_{\mathbb{V}_i, \mathbb{V}_j} \leq \sum_{(\mathbb{V}_i, \mathbb{V}_j) \in \mathcal{L}_{\tau}}{x}^{\tau}_{\mathbb{V}_i, \mathbb{V}_j}, \label{c_y_upper_bound_constraints}
	\end{equation}
	where $M$ is a big constant commonly used in logical constraints reformulation in integer programming, and can be set as $10^6$. It can be checked that if ${x}^{\tau}_{\mathbb{V}_i, \mathbb{V}_j} = 0$, then ${y}^{\tau}_{\mathbb{V}_i, \mathbb{V}_j} =0$; otherwise, ${y}^{\tau}_{\mathbb{V}_i, \mathbb{V}_j} <=M$ automatically holds and has no restriction on ${y}^{\tau}_{\mathbb{V}_i, \mathbb{V}_j}$.
Additionally, for source GT and destination GT, there is
\begin{equation}
\sum_{\mathbb{V}_k \in \mathcal{V}^{\tau} - \{\mathbb{V}_s\}} {y}^{\tau}_{\mathbb{O}_s, \mathbb{V}_k} = 1. \label{c_source_node_y0_constraints}
\end{equation}
\begin{equation}
\sum_{\mathbb{V}_k \in \mathcal{V}^{\tau} - \{\mathbb{O}_d\}} {y}^{\tau}_{\mathbb{V}_k, \mathbb{O}_d} = \sum_{(\mathbb{V}_i, \mathbb{V}_j) \in \mathcal{L}_{\tau}}{x}^{\tau}_{\mathbb{V}_i, \mathbb{V}_j}. \label{c_sink_node_y0_constraints}
\end{equation}

\subsubsection{\text{Sub-tour elimination constraints}}
If the path passes through a forwarding node only once, the order of its outgoing link is the next hop of its incoming link, thus $ \sum_{\mathbb{V}_k:(\mathbb{V}_{\xi},\mathbb{V}_k) \in \mathcal{L}^{\tau}} {y}^{\tau}_{\mathbb{V}_{\xi}, \mathbb{V}_k} = \sum_{\mathbb{V}_k:(\mathbb{V}_k, \mathbb{V}_{\xi}) \in \mathcal{L}^{\tau}} {y}^{\tau}_{\mathbb{O}_k, \mathbb{V}_{\xi}} + 1$.
If the path passes through a forwarding node twice, the sum of the orders of its outgoing links is larger than the sum of the order of its incoming links by $2$, which is $ \sum_{\mathbb{V}_k:(\mathbb{V}_{\xi},\mathbb{V}_k) \in \mathcal{L}^{\tau}} {y}^{\tau}_{\mathbb{V}_{\xi}, \mathbb{V}_k} = \sum_{\mathbb{V}_k:(\mathbb{V}_k, \mathbb{V}_{\xi}) \in \mathcal{L}^{\tau}} {y}^{\tau}_{\mathbb{O}_k, \mathbb{V}_{\xi}} + 2$.
Two different cases can be combined together as one equation, namely, $ \forall \mathbb{V}_{\xi} \in \mathcal{V}^{\tau} -\{\mathbb{O}_s, \mathbb{O}_d\}$,
\begin{equation}
\sum_{\mathbb{V}_k:(\mathbb{V}_k, \mathbb{V}_{\xi}) \in \mathcal{L}^{\tau}}( {y}^{\tau}_{\mathbb{O}_k, \mathbb{V}_{\xi}} + {x}^{\tau}_{\mathbb{O}_k, \mathbb{V}_{\xi}}) = \sum_{\mathbb{V}_k:(\mathbb{V}_{\xi},\mathbb{V}_k) \in \mathcal{L}^{\tau}} {y}^{\tau}_{\mathbb{V}_{\xi}, \mathbb{V}_k}. \label{c_path_sub_tour_elim_constraints}
\end{equation}

	\subsection{Problem formulation}
	The objective is to minimize the end to end delay of the path satisfying the application $\mathcal{A}$'s requirements.
	Therefore, the problem can be formulated as follows:
	\begin{equation*}
	\begin{split}
	\mathbf{P1: } \mathop  \text{min }  & \sum_{(\mathbb{V}_i, \mathbb{V}_j) \in \mathcal{L}^{\tau}} {x}^{\tau}_{\mathbb{V}_i, \mathbb{V}_j} \cdot {D}^{\tau}_{\mathbb{V}_i, \mathbb{V}_j}
	\\  \text{ s.t. }  & (\ref{c_source_node_constraints}) - (\ref{c_path_sub_tour_elim_constraints}).
	\end{split}
	\end{equation*}

	The problem \textbf{P1} is an integer linear programming (ILP) problem, as both the objective and constrains are linear functions. Such problem can be solved by commercial integer programming solvers, such as Gurobi \cite{gurobi}, using the classical branch and bound (B\&B) method. However, for even moderate scale SNs, searching for the optimal solution using the B\&B is still high, because the computational complexity of the B\&B is mainly related to the number of binary variables, the number of total constraints and the scale of the network. In particular, the total number of binary variables in \textbf{P1} is $|\mathcal{L}^{\tau}|$, and the total number of constraints is $H = 4 + (|\mathcal{V}^{\tau}|-2) \cdot 3 + 2 + |\mathcal{L}^{\tau}| = 3 |\mathcal{V}^{\tau}| + |\mathcal{L}^{\tau}| $. According to \cite{wolsey1999integer}, the B\&B method for solving \textbf{P1} has the worst case time complexity of $\mathcal{O}(|\mathcal{L}^{\tau}| \cdot H^{|\mathcal{L}^{\tau}|})$, which is exponential with the total number of links $|\mathcal{L}^{\tau}|$. Even for a medium scale \textbf{P1}, the running time can be in the order of minutes to hours, or even days for larger network sizes.
	Therefore, it is necessary to exploit the special structure of \textbf{P1} and obtain more efficient methods.

	\section{The proposed graph-based algorithms}
Basically, solving \textbf{P1} is equal to finding the minimum delay $(\mathbb{O}_s-\mathbb{O}_d)$ path with capacity no less than $C_a$ and including at least one satellite deployed with function $f_a$.
In the following, firstly, we remove some capacity insufficient links and adopt the k-shortest path (KSP) algorithm\textcolor{black}{\cite{yen1971finding}} to solve it. After analyzing its drawbacks, an alternative algorithm is proposed by exploiting the special structures of the solution space.

\subsection{KSP-based VF constrained simple path algorithm}

The main idea is to iteratively find the $k$-th shortest delay paths from source GT $\mathbb{O}_s$ to destination GT $\mathbb{O}_d$ in the residual snapshot graph removing capacity-insufficient links, and stop until the required path is found or the $k$-th path violates the required path delay bound $D_a$. Firstly, all the links with capacity less than $C_a$ are removed to trim the solution space.
Then, the searching procedure starts from the shortest delay path and each iteration a $k$-th shortest path $p_k$ is generated and checked.  If the new generated path $p_k$ satisfies the delay bound requirement, it will be checked whether it contains a satellite supporting function $f_a$. If a feasible path is found, the searching procedure terminates. Otherwise, it will stop until the delay of the $k$-th path exceeds the required delay $D_a$.
The detailed KSP-based algorithm is shown in Algorithm 1.

\begin{algorithm}
	\caption{KSP-based VF constrained simple path algorithm}

	\label{array-sum}
	\hspace*{0.02in} {\bf Input:}

	$ \mathcal{G}^{\tau}=\{\mathcal{V}^{\tau}, \mathcal{L}^{\tau} \}$,
	and $\mathcal{A} = \{\mathbb{O}_s, \mathbb{O}_d, f_a, C_a, D_a \}$.  \\
	\hspace*{0.02in} {\bf Output:}
	The virtual function constrained simple path $p_{*}$.
	\begin{algorithmic}[1]
		\item  Initialize $\mathcal{G}^{\tau}_{C_a} \leftarrow \mathcal{G}^{\tau}$.
		\item  \textbf{for} each link $(\mathbb{V}_i, \mathbb{V}_j) \in \mathcal{L}^{\tau}$ \textbf{do}
		\item  \ \ \ \  \textbf{if} $C^{\tau}_{(\mathbb{V}_i, \mathbb{V}_j)} < C_a $ \textbf{then}

		\item  \ \ \ \ \ \ \ \   \text{Remove} link $(\mathbb{V}_i, \mathbb{V}_j) $ from $\mathcal{G}^{\tau}_{C_a} $.
		\item $k=0$, $p_{*} = \varnothing $, $D_{p_0} = 0 $.

		\item \textbf{while} $D_{p_k} \leq D_a $ \textbf{do}
		\item 	\hspace*{0.02in} \hspace*{0.02in}
		$k \longleftarrow k+1$
		\item 	\hspace*{0.02in} \hspace*{0.02in}
		Find the $k$-th shortest path $p_k$ from $\mathbb{O}_s$ to $\mathbb{O}_d$ in $\mathcal{G}^{\tau}_{C_a}$.
		\item 	\hspace*{0.02in} \hspace*{0.02in}
		Calculate the delay of the $k$-th path $p_k$ as $D_{p_k}$.
		\item  \hspace*{0.02in} \hspace*{0.02in} \textbf{if} $\exists \mathbb{S}_{\xi} \in p_k, w^{f_a}_{\mathbb{S}_{\xi}}\geq 1 $ \textbf{then}
		\item  \hspace*{0.02in} \hspace*{0.02in} 	\hspace*{0.02in} \hspace*{0.02in} $p_{*} = p_k$, break.
		\item \textbf{return} $p_{*}$.
	\end{algorithmic}
\end{algorithm}

Obviously, the running time of \text{Algorithm 1} depends on how many paths within delay bound $D_a$ need be iterated before a path containing a satellite deployed with function $f_a$ is found.
The drawbacks of Algorithm 1 are summarized as follows:

\begin{itemize}

	\item  \textit{Sub-optimality:} Since the k-shortest path algorithm computes simple paths with no repeating nodes, the optimal path which have repeated satellites can not be obtained by this algorithm. Even worse, such a drawback may cause applications to be rejected though there is more than one feasible non-simple path with repeated nodes.

	\item \textit{Instability:} When the delay bound $D_a$ is small, the algorithm will stop quickly since there are fewer paths to be checked. However, when $D_a$ is large, it may take a very long time for the algorithm to iterate paths. Therefore, the running time of the algorithm is related to both the value of $D_a$ and the number of function enabled satellites, which makes the algorithm unstable and pseudo-polynomial.

	\item  \textit{Non Scalability: } Although a single path can be found in polynomial time, e.g., $\mathcal{O}(|\mathcal{L}^{\tau}| + |\mathcal{V}^{\tau}|\log|\mathcal{V}^{\tau}|) $ by using the Dijkstra's algorithm, the number of infeasible paths in a densely connected graph with thousands of nodes can be extremely large, i.e., $\mathcal{O}(|\mathcal{V}|!)$ in a complete graph of order $|\mathcal{V}^{\tau}|$. For large scale SNs, the KSP-based method will be very time consuming.
\end{itemize}

\textit{Complexity analysis:}
In fact, based on the KSP algorithm proposed by \cite{yen1971finding}, finding $K$ number of simple paths requires $O(K|\mathcal{V}^{\tau}|^3)$ operations. However, in the worst case, given a complete graph of order $|\mathcal{V}^{\tau}|$ with no function enabled satellite and $D_a = +\infty$, the computational complexity of Algorithm 1 is $O(|\mathcal{V}^{\tau}|! \cdot |\mathcal{V}^{\tau}|^3)$, which is prohibitively time consuming and can not be deployed in large scale SNs.

\subsection{The proposed VF-aware shortest path algorithm}

	In practice, both the ILP-based method and KSP-based method have drawbacks, preventing these two solvers from being applicable in large-scale SNs.
	To address this bottleneck, we propose an alternative method by bidirectional path seeking from source and destination to the functional satellites (satellites can providing required functions), such that the problem can be solved optimally by running two times of Dijkstra's algorithm combined with a node selection procedure.

	Instead of directly searching for candidate paths, we collect the functional satellites during the link filtering operation (same as Algorithm 1) before calculating the path.
	As a result, the path seeking problem of \textbf{P1} can be decomposed into two sub-problems: 1) Seeking the shortest path from source GT to every functional satellites; 2) Seeking the shortest path from every functional satellites to destination GT.
	Since there can be multiple functional satellites, the complexity of path calculation is too high.
	Fortunately, the Dijkstra's algorithm has the special property that the computation of one single-source single-sink shortest path can obtain all the shortest paths from the single-source to all other destinations. In other words, the complexity of single-source single-sink shortest path is the same as that of the single-source multi-sink shortest path.
	Therefore, we view the functional satellites as multiple sinks, and run two times of Dijkstra's algorithm, one in the forward direction from source GT to functional satellites, and another in the reverse direction from the destination GT to all functional satellites by reversing all links of the original snapshot graph. Finally, the VF-aware shortest path can be found by joining the partial paths, which is the optimal solution of \textbf{P1}.

\begin{algorithm}
	\caption{The VF-aware shortest path (VFSP) algorithm}
	\label{array-sum}
	\hspace*{0.02in} {\bf Input:}
	$\mathcal{G}^{\tau}=\{\mathcal{V}^{\tau}, \mathcal{L}^{\tau} \}$, and $\mathcal{A} = \{\mathbb{O}_s, \mathbb{O}_d, f_a, C_a, D_a \}$.\\
	\hspace*{0.02in} {\bf Output:} The virtual function constrained shortest path $p^*$.

	\begin{algorithmic}[1]
		\item  Initialize $\mathcal{G}^{\tau}_{C_a} = \varnothing$, $\mathcal{V}^{\tau}_{f_a} = \emptyset$, $p_{*} = \varnothing$, and $D_{p_{*}} = D_a + 1 $.
		\item  \textbf{for} each link $(\mathbb{V}_i, \mathbb{V}_j) \in \mathcal{L}^{\tau}$ \textbf{do}
		\item  \ \ \ \  \textbf{if} $C^{\tau}_{(\mathbb{V}_i, \mathbb{V}_j)} \geq  C_a $ \textbf{then}

		\item  \ \ \ \ \ \ \ \   \text{Add} link $(\mathbb{V}_i, \mathbb{V}_j) $ to $\mathcal{G}^{\tau}_{C_a} $.
		\item  \ \ \ \ \ \ \ \  \textbf{if} $w^{f_a}_{\mathbb{V}_{i}}\geq 1 $ \textbf{then}

		\item  \ \ \ \ \ \ \ \  \ \ \ \   $\mathcal{V}_{f_a} \leftarrow \mathcal{V}_{f_a} \cup \{\mathbb{V}_i\}$.
		\item  \ \ \ \  \ \ \ \  \textbf{if} $w^{f_a}_{\mathbb{V}_{j}}\geq 1 $ \textbf{then}

		\item  \ \ \ \ \ \ \ \ \ \ \ \   $\mathcal{V}_{f_a} \leftarrow \mathcal{V}_{f_a} \cup \{\mathbb{V}_j\}$.
		\item  \textbf{if} $\mathcal{V}_{f_a} = \emptyset $ \textbf{then}
		\item  \ \ \ \  \textbf{return} -1. // No feasible path.
		\item  Reverse all links in $\mathcal{G}^{\tau}_{C_a}$ to build a reversed graph $\bar{\mathcal{G}}^{\tau}_{C_a}$.
		\item  Calculate shortest paths from  $\mathbb{O}_s$ to all other nodes in $\mathcal{G}^{\tau}_{C_a}$.
		\item  Calculate shortest paths from  $\mathbb{O}_d$ to all other nodes in $\bar{\mathcal{G}}^{\tau}_{C_a}$.
		\item  \textbf{for} each satellite $\mathbb{S}_i \in \mathcal{V}^{\tau}_{f_a} $ \textbf{do}
		\item  \ \ \ \  \textbf{if} $D_{p_{(\mathcal{O}_s, \mathbb{S}_i)}} + D_{p_{(\mathcal{O}_d, \mathbb{S}_i)}} \leq D_{p_{*}}$ \textbf{then}

		\item  \ \ \ \ \ \ \ \   $D_{p_{*}} = D_{p_{(\mathcal{O}_s, \mathbb{S}_i)}} + D_{p_{(\mathcal{O}_d, \mathbb{S}_i)}}$, $p_{*} =p_{(\mathcal{O}_s, \mathbb{S}_i)} + \bar{p}_{(\mathcal{O}_d, \mathbb{S}_i)} $.
		\item  {\bf return}  $p_*$.
	\end{algorithmic}
\end{algorithm}
\textit{Complexity analysis:}
Based on the VFSP, finding all shortest paths and functional satellites requires $O(|\mathcal{L}^{\tau}|)$ operations. And the computational complexity for Dijkstra's algorithm for computing one source multiple sink shortest path is $\mathcal{O}(|\mathcal{L}^{\tau}|+|\mathcal{V}^{\tau}| \log{|\mathcal{V}^{\tau}|})$. Therefore, the worst case computational complexity of Algorithm 2 is $\mathcal{O}(3|\mathcal{L}^{\tau}| + (2\log{|\mathcal{V}^{\tau}|}+1) |\mathcal{V}^{\tau}|)$, which is polynomial with the increase of network size.

\section{Evaluation}
\subsection{Simulation setup}
We conduct the simulation based on the starlink constellation, which is the largest LEO satellite system by far.
 Specifically, we randomly choose $100$-$2300$ satellites from 2694 active starlink satellites in standard object database of systems tool kit (STK),
a third-party software that updates continuously to simulate the movement of real-world satellites.
The ground terminals are randomly distributed and located in around Xi'an ( $34.27^{\circ} N, 108.93^{\circ} E$ ), Beijing ($40^{\circ} N, 116^{\circ} E$), Sanya ( $18^{\circ} N, 109.5^{\circ} E$) and Kashi ($39.5^{\circ} N, 76^{\circ} E$).
The contact plans including communication opportunities of each pair of satellites in the studied network are obtained by using \textit{compute access} function in STK.
The transmission rate of both ISLs and GSLs is uniformly selected from $[300,350]$ Mbps as in \cite{fu2020remote}.
The link delays of both ISLs and GSLs are in the range of $[5, 15]$ ms as in \cite{chen2022robust}.

The simulation time horizon $T$ is from \text{2022-8-9 04:00:00} to \text{2022-8-10 04:00:00}, and is split into multiple time windows by the time division mechanism in \cite{wang2022enhancing}.
All the three algorithms are implemented using Python programming language.

\subsection{Simulation results and analysis}
\subsubsection{Running time performance with increasing satellite network sizes}
Fix the deployed function number $|\mathcal{F}|=3$ and the percentage of function enabled satellites as $10\%$, we consider 5000 random applications. Each application involves a pair of GTs, and requires the transmission capacity and the end to end delay randomly selected from $[5, 100]$ Mbps, and $[20, 150]$ ms, respectively.
\textcolor{black}{Figure 2} reports the running time of three different algorithms with the number of satellites varying from $100$ to $2300$, with a step-wise increase of $100$.
As the network size increases, the average running time of all the three algorithms increases.
However, both the graph-based methods (KSP, VFSP) are significantly faster than the ILP method.
This can be explained as follows. For the ILP-based method, as the satellite number becomes large, more transmission edges are added into the snapshot graph, resulting in increased number of decision variables and exponential growing running times. Similarly, the increased transmission opportunities also introduce more infeasible paths in the snapshot graph for the KSP-based method to explore, thus the number of iterations increases. Since the computation complexity of the proposed VFSP mainly depends on the function enabled satellite numbers in a network, its running time increases linearly with the increase of satellite numbers.

\begin{figure}[t!]
	\centering
		\includegraphics[width=60mm]{Graph6-eps-converted-to.pdf}\\
	\caption{Running times versus different satellite numbers. }
\end{figure}

\begin{figure}[t!]
	\centering
			\includegraphics[width=60mm]{Graph13-eps-converted-to.pdf}\\
	\caption{Acceptance ratios versus function enabled satellite numbers. }
\end{figure}

\subsubsection{Network performance with increasing function enabled satellites}
Fix the satellite number as $1000$, we vary the percentage of function enabled satellites to investigate the impact of the function deployment ratio on network performance. In Figure 3, the acceptance ratio versus the percentage of function enabled satellites is plotted from $5\%-100\%$ with an increase step of $5\%$, where $5\%$ means that $50$ out of $1000$ satellites are deployed with virtual functions. For comparison, each function enabled satellite accommodates at most one application. As seen from Figure 3, more function enabled satellites can bring higher acceptance ratio for both the KSP and VFSP. It is as expected that the acceptance ratio of the VFSP is higher than that of the KSP, since the VFSP can explore more non-simple paths to accommodate applications while the KSP can not. Moreover, when the percentage of function enabled satellites is less than $20\%$, the acceptance ratio grows faster since the function resources are the main bottleneck of the network.

Figure 4 and Figure 5 plot the average delay and average number of hops of the paths with the increase of the percentage of function enabled satellites, respectively. The path delay of the KSP is smaller than that of the VFSP, which is expected since applications with no feasible simple paths are rejected by the KSP.
However, such applications can be fulfilled by the VFSP with non-simple paths with repeated nodes thus longer delays. The trend in Figure 5 is similar to that of Figure 4, as smaller end to end delays correspond to fewer hop numbers. However, the lower delay and fewer hops performance of the KSP method is at the cost of reduced acceptance ratios, since a large proportion of feasible applications can be rejected.

\begin{figure}[t!]
	\centering
	\includegraphics[width=60mm]{Graph11-eps-converted-to.pdf}\\
	\caption{Path delays versus function enabled satellites. }
\end{figure}

\begin{figure}[t!]
	\centering
	\includegraphics[width=60mm]{Graph12-eps-converted-to.pdf}\\
	\caption{Number of hops versus different satellite numbers. }
\end{figure}

	\section{Conclusion}

In this work, we investigate the routing strategy for real time applications in large-scale SNs with VFs. We identify that such a routing problem can be formulated as an integer linear programming problem, which incurs exponential complexity by using the branch and bound method. To solve the problem efficiently, two alternative graph-based algorithms from different standpoints are proposed, analyzed and compared, where the latter one can obtain the optimal solution in polynomial time with a low computation complexity and more stable performance in large scale networks. Simulations conducted on starlink constellation with thousands of satellites verify the performance of the proposed algorithms.

	%

	\ifCLASSOPTIONcaptionsoff
	\newpage
	\fi

	%

	\bibliographystyle{IEEEtran}