\section{Introduction }
Topological materials are among the most intensely researched current topics, and encompass not just topological insulators~\cite{kane2005quantum,kane2005z,PhysRevLett.96.106802,bernevig2006quantum, fu2007topological,PhysRevLett.106.106802, hsieh2012topological,tanaka2012experimental,dziawa2012topological,PhysRevB.78.195125}
but also semi-metals with topological nontrivial nodal structures~\cite{burkov2016topological,doi:10.1146/annurev-matsci-070218-010049, doi:10.7566/JPSJ.87.041001, fang2016topological,doi:10.1146/annurev-matsci-070218-010023}. Such nodal materials have provoked widespread studies due to their unconventional density of states and bandstructure, which have lead to new or enhanced avenues for non-linear optical responses, electron tunneling behavior, higher harmonic generation, superconductivity and quantum Hall effects~\cite{RevModPhys.82.3045,doi:10.1146/annurev-conmatphys-031214-014501,RevModPhys.91.015006, ozawa2019topological, doi:10.7566/JPSJ.87.041001,dong2010quantum,pang2015evidence,wang20173d,molina2018surface,zhang2019cyclotron,zhang2019quantum,ang2020universal,bhattacharyya2020evidence}. The potentially intricate topology of nodal loops in 3D space have also inspired their design and realization in metamaterials and lossy or non-reciprocal media~\cite{lu2013weyl,xiao2015synthetic,lin2016photonic,yang2016acoustic,chen2016photonic,chang2017multiple,noh2017experimental,lin2017line,li2018weyl,yang2018ideal,gao2018experimental,lee2018tidal,xia2019observation,zhou2019exceptional,lee2020imaging,cerjan2019experimental,li2019emergence,jia2019observation}.

Despite their seemingly exotic nature, various Nodal line semimetal (NLSM) materials have recently been experimentally characterized, such as PbTaSe$_2$~\cite{bian2016topological, PhysRevB.93.245130}, BiTeI~\cite{crepaldi2012giant}, Mg$_3$Bi$_2$~\cite{doi:10.1021/acs.jpclett.7b02129,chang2019realization,zhou2019experimental}, ZrSiTe and ZrSiSe~\cite{hu2016evidence,PhysRevX.10.011026}, ZrSiS~\cite{schoop2016dirac,Fueaau6459,PhysRevB.100.085137}, BaAgAs~\cite{mardanya2019prediction}, TaN~\cite{PhysRevB.93.241202},  Ca$_3$P$_2$~\cite{chan20163,xie2015potential},  SrAs$_3$~\cite{song2020photoemission,li2018evidence}, CaAgX (X=P,As)~\cite{yamakage2016line,takane2018observation}, ZrGeX$_\text{c}$ (X$_\text{c}$= S, Se, Te) family~\cite{nakamura2019evidence, guo2019electronic}, magnetic semimetals EuB$_6$~\cite{PhysRevLett.124.076403}, spin gapless semimetals~\cite{zhang2020nodal}, the centrosymmetric superconductor SnTaS$_2$ ~\cite{chen2019superconducting} and  PbTaS$_2$~\cite{gao2020superconducting}. In the form of closed loops, they have also been observed as nodal links in CaAuAs~\cite{PhysRevB.98.085122}, nodal chains in TiB$_2$ ~\cite{yi2018observation, liu2018experimental} and nodal line networks in RuO$_2$~\cite{jovic2019dirac}. Other materials such as Ti$_3$X (X=Al, Ga, Sn, Pb) family~\cite{PhysRevB.97.235150}, YH$_3$\cite{shao2018nonsymmorphic}, YoCoC$_2$~\cite{xu2019topological}, MnN~\cite{PhysRevMaterials.3.084201}, TiTaSe$_2$~\cite{bian2016drumhead}, CaP$_3$ family~\cite{PhysRevB.95.045136}, ABC-stacked graphdyiyne~\cite{graphdiyne}, layered X$_2$Y (X=Ca, Sr, Ba; Y=As, Sb, Bi)~\cite{niu2017two} have also been theorically proposed to possess
nodal loops.

In some of these NLSMs, the conduction and the valence bands intersect to form potentially very intricate structures, even more intricate than common examples in knot theory and prototypical nodal knot setups~\cite{AlexanderPolynomial,vortexKnot,dennis2017constructing,lee2020imaging,bode2019constructing}. For instance, the material Co$_2$MnGa was theorized~\cite{CMG_theory} to possess a complicated network of 3D band crossings characterized by various types of non-trivial nodal linkages and coupled chains enabled by perpendicular mirror planes. This was subsequently confirmed experimentally~\cite{belopolski2019discovery}. Another proposal of interlocking nodal chains may be realized in carbon networks, consisting of armchair graphene nanoribbon~\cite{LI2020563}. While transport and optical response properties are already well-studied for simpler NLSMs ~\cite{doi:10.1146/annurev-matsci-070218-010023,PhysRevLett.124.076403,PhysRevB.100.085137,yang2019conductivity,guo2019electronic, PhysRevB.99.241102,PhysRevB.95.245113,PhysRevB.102.035164}, those of more complicated nodal \emph{loop} semimetals are still not well understood. Recently, it was proposed that topological nodal linkages can significantly enhance optical response non-linearity and hence HHG (higher harmonics generation)~\cite{PhysRevB.102.035138}. However, the non-linear response does not afford any topological quantization, unlike linear response via the Kubo formula, and a complete understanding of the response properties of realistic nodal loop material necessitates a systematic study of how the geometry of the nodal structure, together with its topology, interplay and lead to various anisotropic and non-linear behavior.

As such, this work shall be concerned with a systematic investigation on how the relative shapes and linkage of nodal loops can lead to various anisotropic and non-linear responses components, and how the various responses come together into a bigger picture that reveals the overall nodal topology and geometry. Following a review of semi-classical response theory in Sect. II, we introduce a canonical model of two nodal loops whose linkage and shapes can be independently tuned. This model serves as an abstraction of the multiple nodal touchings and linkages in realistic nodal materials, for instance Co$_2$MnGa. In Sect. III, we study the various components in the response tensor of our model system, some whose non-linear properties have never been investigated. In Sect. IV, we show how these results help piece together a response surface whose evolution with field strength encapsulates the full information about the response anisotropy and non-linearity, and to some extent the nodal structure and its dispersion. For our model, the direction which the response surface points towards in the high field limit depends on whether the nodal structure is topologically linked.

\section{Non-linear Semi-Classical Optical Response}
To systematically relate the Hopf link response to its detailed momentum-space profile, we first review the theory of non-linear semi-classical response, and next introduce our model with tunable nodal linkage. Similar semi-classical approaches have been highly successful in explaining phenomena like Hall effects and quantum oscillations in diverse settings~\cite{xiao2010berry,price2016measurement,petrides2018six,lee2018electromagnetic,lohse2018exploring}, as well as Bloch oscillations and Berry curvature effects in the context of HHG~\cite{schubert2014sub,liu2017high}. HHG refers to the high harmonics generation induced by strong non-linear interactions when a very intense laser pulse is focused into a sample.

\subsection{Semi-Classical Response Theory}
The non-linear response of nodal materials emerges even at the semi-classical level, where many-body processes are simply encapsulated by a non-equilibrium occupation function $F(\varepsilon(\mathbf{k}-e\mathbf{A(t)}))$ that weighs different contributions to the response optical current~\cite{mikhailov2007non}:
\begin{equation}
\mathbf{J}=\int F(\varepsilon(\mathbf{k}-e\mathbf{A(t)}))\langle\mathbf{k}|\mathbf{\hat J}|\mathbf{k}\rangle d^3\mathbf{k}.
\label{J}
\end{equation}
Here, $F(\varepsilon)=(1+e^{\beta(\varepsilon(\mathbf{k})-\mu)})^{-1}$ is the \emph{equilibrium} Fermi-Dirac occupation function that depends on the energy dispersion $\varepsilon(\mathbf{k})$, and $e\mathbf{A}(t)=e\int_{-\infty}^t\mathbf{E}(t')dt'$ is the impulse on an electron $e$ due to an external electric field $\mathbf{E}(t)$. In frequency space, we can always replace $\bold A$ with $\bold E/i\Omega$. The current operator for a given Hamiltonian $\hat H$ is $\mathbf{\hat J}=\frac{\partial\hat H}{\partial\mathbf{k}}$, which reduces to the dispersion velocity $\langle\mathbf{k}|\mathbf{\hat J}|\mathbf{k}\rangle=\bold v(\mathbf{k})=\frac{\partial \varepsilon(\mathbf k)}{\partial\mathbf{k}}$ in the translation-invariant intra-band case. We express this optical response in material-dependent units of $nev_F\sim10^{12}$ Am$^{-2}$ where $n$ and $v_F$ are the charge carrier density and the Fermi velocity. In real materials, we have $v_F$ to be $2.22\times10^5$ ms$^{-1}$ and $4.0\times10^5$ ms$^{-1}$ in ZrGeSe ~\cite{doi:10.1063/1.5084090} and CaAgAs~\cite{PhysRevB.95.245113} respectively, $n$ to be $3.37\times10^{26}$ m$^{-3}$, $2.98\times10^{26}$ m$^{-3}$ and $1.7\times10^{26}$ m$^{-3}$ in YbCdGe~\cite{PhysRevB.99.241102}, CaCdSn~\cite{PhysRevB.102.035164} and CaAgAs~\cite{PhysRevB.95.245113} respectively. In these NLSMs, the Fermi surface is a simple nodal loop. The optical response of real materials with more sophisticated nodal structure, beyond the simplest nodal loop, have also been computed in Ref. \cite{PhysRevB.102.035138}.

Eq.~\ref{J} holds in the ballistic limit, where it is the exact solution to the semi-classical Boltzmann equation and the semi-classical equations of motions for electronic wavepackets~\cite{mikhailov2007non}. This requires $\Omega\tau\gg 1$, where $\Omega$ is the optical frequency defined by $\mathbf{E}(t)\sim\mathbf{E}(0)e^{i\Omega t}$, and $\tau$ is the relaxation time due to impurity scattering.  For a nodal material with relaxation time $\tau\sim 10^{-12}$ to $10^{-13}\,s$, comparable to that in high quality Graphene samples~\cite{TransportTimeGraphene}, $\Omega\tau\gg 1$ can be achieved in the teraHertz regime of $\Omega \sim 50$-$100\,$THz. The relaxation times can be computed from a microscopic model for the scattering processes~\cite{ho2018theoretical}. In this ballistic regime, scattering processes cannot catch up with the much shorter oscillation timescales~\cite{liu2017high,ndabashimiye2016solid}.

According to Eq.~\ref{J}, the response current arises from contributions at $\bold k$ for which the Fermi-Dirac occupation $F(\varepsilon(\mathbf{k}))$ is non-vanishing. As such, it strongly depends on the shape of the (occupied) Fermi region is the Brillouin zone (BZ), particularly its co-dimensionality. This intra-band response is the expectation of $\bold v(\bold k+\bold e \bold A(t))$ within the occupied region that is dynamically shifted by the external electromagnetic impulse. When the occupied region is not a ``blob'' in the BZ but a thin Fermi ``tube'' of nontrivial co-dimension, the expectation of $\bold v(\bold k+\bold e \bold A(t))$ can fluctuate wildly with $\bold A(t)$ due to ``destructive interference'' of $\bold v$ at different $\bold A(t)$ shifts, as previously pointed out in Refs.~\cite{mikhailov2007non,lee2015negative,PhysRevB.102.035138}. In Ref.~\cite{mikhailov2007non}, this observation was first used to explain the non-linear response of Graphene due to its vanishing density of states at small chemical potential. This was generalized in~\cite{lee2015negative} to nodal loop structures i.e. where two bands intersect along a loop in the BZ, where a stronger non-linearity from a characteristic non-linear response curve was reported. In Ref.~\cite{PhysRevB.102.035138}, this non-linear response was further shown to be strongly enhanced in a 3D material with non-trivial nodal topology i.e. when band intersections form loops that are topologically linked.

To more concretely understand why nodal structures with nontrivial linkages have strong non-linear responses, consider small chemical potentials $\mu$ such that the occupied states take the form of nodal ``tubes'' along nodal lines, which are then shifted by the electromagnetic impulse $\int^t\mathbf{E}(t')dt'$ away from the original nodal structure in the BZ. By differentiating Eq.~\ref{J} with respect to the vector potential, we obtain~\cite{PhysRevB.102.035138} for small $\mu$,
\begin{equation}
\frac{\partial J_i}{\partial A_j}\approx2\sum_{\alpha\in\text{NLs}}\frac{\mu}{\hat{e_j}\cdot\mathbf{v_F}}\frac{\partial^2\varepsilon(\mathbf{k}^\alpha+e\mathbf{A})}{\partial k_i \partial k_j}
\label{dJdA}
\end{equation}
where $\mathbf{k}^\alpha$ labels the trajectories of all the nodal lines $\alpha$ in the BZ and $\hat{e}_j\cdot\mathbf{v}_F$ is the component of the Fermi velocity of the $\alpha$-th NL at momentum $\mathbf{k}^\alpha$ along the applied field. In Eq.~\ref{dJdA}, the differential (optical) response tensor of a nodal line structure is given by a sum over the Hessian of the dispersion at the nodal lines shifted by $\mathbf{A}$ in the BZ. Evidently, we expect the  response to be non-linear whenever $\frac{\partial J_i}{\partial A_j}$ is significantly non-constant, which has to be the case around nodal crossings with vanishing gap. In particular, in the diagonal sector where $i=j$, $\varepsilon(\mathbf{k}^\alpha+e\mathbf{A})$ has to pass through a singularity whenever two nodal loops are linked in the direction $\bold A$, since the locus of $\mathbf{k}^\alpha+e\mathbf{A}$ for one of the loops $\alpha$ must intersect another loop $\alpha'$ as $\bold A$ increases. This leads to the enhancement of diagonal response non-linearity in topologically linked nodal loops, which can be quantified by the extent of HHG it causes.

However, the abovementioned argument for non-linearity enhancement pertains only to the diagonal $i=j$ sector, where the response is probed in the direction where the applied field displaces the nodal loops. It remains an open question if a nodal linkage also possess signature response behaviors in the other directions, in both longitudinal and transverse sectors. In this work, we shall focus on elucidating the signature contributions to the optical response due to a generic nodal linkage, particularly on how non-topological contributions from the dispersion profile competes with those from the topological linkage. While real nodal materials may also experience response contributions from the Berry curvature and inter-band scattering, these additional contributions are very material-dependent, and do not form the focus of this work, which aims to explore the most generic response behavior expected from nodal linkages.

\subsection{Canonical Hopf link Model}
In this work, we have chosen to focus on the simplest topologically non-trivial nodal linkage, which is the Hopf link - a nodal configuration consisting of two nodal loops that will be linked when they are sufficiently close to each other. By introducing a canonical two-band Hopf link model where the key parameters $r_y$, $r_x$ and $d$ (Fig.~\ref{Hopf_shape}a) are all independently adjustable, we can isolate the various geometric and topological factors that influence its non-linear response. The parameters are:
\begin{itemize}
    \item $r_y$: width of the loops in the $k_y$-direction or longitudinal direction along which the two loops are separated;
    \item $r_x$: width of the loops in the $k_x$ or $k_z$ transverse directions;
		\item $2d$: the displacement between the two loops, centered at $(0,\pm d,0)$. They are linked if $r_y>d$, since the longitudinal separation between them is given by $2(r_y-d)$ (Fig.~\ref{Hopf_shape}).
\end{itemize}

\begin{figure}
\centering
\includegraphics{Fig1.pdf}
    \caption{(a) The Fermi surface of the Hopf link in $\mathbf{k}$-space, with the parameters $r_x$, $r_y$, $r_z$ and $d$ as labeled. We assumed $r_x=r_z$. The quantity $2(r_y-d)$ controls the extent of linkage in the Hopf link. (b) The three distinct possibilities of the Hopf structure are: Hopf link ($r_y<d$), nodal chain ($r_y\approx d$) and unlinked nodal loops ($r_y>d$), illustrated here for a non-vanishing $\mu=0.1$.}
    \label{Hopf_shape}
\end{figure}

Our two-band nodal Hopf link Hamiltonian takes the form
\begin{eqnarray}H_{\text{hopf}}(\mathbf{k})=h(k_x,k_y,k_z)\sigma_x+g(k_x,k_y,k_z)\sigma_y
\end{eqnarray}
where the functions $h(\mathbf{k})$ and $g(\mathbf{k})$ are obtained by deforming a well-known Hopf link model derived through the Hopf map~\cite{PhysRevB.96.041202,PhysRevB.97.155140,PhysRevLett.118.147003,lee2020imaging}, for which the Hopf link parameters cannot be easily adjusted independently~\footnote{An analogous approach~\cite{lee2017band} with complex singularities also does not allow for independent tuning of parameters.}. The gap $2\epsilon(\mathbf{k})=2\sqrt{g(\mathbf{k})^2+h(\mathbf{k})^2}$ closes along two loops given by the simultaneous solutions of $g(\bold k)=0$ and $h(\bold k)=0$. Explicitly, with the intermediate parameters $\xi=\frac{\cos r_y-1}{\cos r_x-1}\frac{1}{\sin d}-1$ and $\Gamma=\frac{\cos r_x\cos r_y -1}{\cos r_x -1}\frac{1}{\sin d }$, $h(\bold k)$ and $g(\bold k)$ are given by
	\begin{widetext}
	\begin{subequations}
	\begin{align}
h(k_x,k_y,k_z)&=\sin^2(k_z-k_x)-(\xi\cos(k_x+k_z)+\cot d \cos(k_y)+\cos(k_z-k_x)-\Gamma)^2-\sin^2(k_x+k_z)+\sin^2(k_y),\\
g(k_x,k_y,k_z)&=2\sin(k_z-k_x)(\xi\cos(k_x+k_z)+\cot d \cos(k_y)+\cos(k_z-k_x)-\Gamma)-2\sin(k_x+k_z)\sin(k_y).
	\end{align}
	\end{subequations}
	\end{widetext}

At chemical potential $\mu=0$, the two nodal loops lie in the planes $k_x=0$ and $k_z=0$, with a common mirror symmetry axis along the $k_y$ line. The nodal loops each touch this symmetry axis at two points, giving rise to a total of four touching points $\pm(r_y\mp d)$. From Fig.~\ref{Hopf_shape}, $r_y-d$ evidently describes the extent of their linkage.  When $r_y>d$, the two nodal loops are linked together, and we call this resulting arrangement to be a Hopf link.

In realistic materials, the chemical potential $\mu$ can often be tuned away from the nodal energy, so as to obtain a finite density of states available for transport. For generic linearly dispersive nodal systems, $\mu$ scales with the thickness of the nodal tube. A large $\mu$ tends to result in thick tubes that inevitably intersects, resulting in a nodal chain even when $r_y\neq d$.

In total, there are three possibilities as illustrated in Fig.~\ref{Hopf_shape}b:
\begin{itemize}
\item Hopf link: The two nodal loops are topologically linked ($r_y>d$ for $\mu\ll 1$),
\item Unlinked nodal loops: The two nodal loops are disjoint and not topologically linked ($r_y<d$ for $\mu\ll 1$),
\item Nodal chain: The two nodal loops touch each other at a point, and are not topologically linked ($r_y\approx d$ for $\mu\ll 1$).
\end{itemize}

Due to the periodicity of its various terms, the Hamiltonian can sometimes admit additional solutions i.e. ``\emph{periodic images}'' within the first BZ. These solutions, analogous to the degenerate valleys in Graphene, occur as two distinct types, as further elaborated in the Appendix. To prevent the confounding ambiguities resulting from the interference of multiple nodal linkages, we shall exclusively study only cases where such periodic images do not exist.

For the rest of this work, we shall set the temperature to a representative value of $10$ K i.e. $\beta:=\frac{eV}{k_BT}=1160$, which can be adjusted to fit the physical temperature of an actual physical scenario by trivially rescaling our canonical model. For definiteness, we choose a chemical potential of $\mu=0.1$ eV, such that the nodal structure consist of thin closed tubes that exhibit distinct crossovers in topology (linked, intersecting, unlinked) as they are shifted due to photon excitations. However, our results shall continue to hold even when typical photon energies are much smaller than $\mu$. We shall study the behavior of the optical response as the three independent parameters $\{d,r_x,r_y\}$ are varied, insofar as the parameters do not admit periodic images, for all three types of cases (linked, unlinked and nodal chain).

\section{Optical Response in various directions}
We next investigate the various diagonal and off-diagonal (Hall) responses of the nodal Hopf link in the principal directions - both longitudinal (along the direction of loop separation $\hat y$) and transverse directions (perpendicular to the direction of loop separation $\hat x$ and $\hat z$).

Consider a sinusoidal time-varying applied electric field signal
$\mathbf{E}(t)=\mathbf{E_0}\cos\Omega t$, which corresponds to the vector potential $\mathbf{A}(t)=\mathbf{A}\sin\Omega t$, where $\mathbf{A}=\mathbf{E_0}/i\Omega$. In the presence of an oscillatory electric field $\mathbf{E}(t)$, the minimally coupled Fermi-Dirac occupation in Eq.~1 will be $F(\varepsilon(\mathbf{k}-e\mathbf{A}(t)))$ which is equivalent to a translation of the Fermi tube of occupied states in the field direction according to the momentum shift $\mathbf{k}\rightarrow\mathbf{k}-e\mathbf{A}(t).$ This picture motivates the $\mathbf{A}$ dependence of $\mathbf{J}$ rather than $\mathbf{E}$. The ballistic criterion $\Omega\sim100$ THz is independent of the nodal energy since the Fermi surface is translationally invariant. Although there are nine possible response components $J_i(A\hat{j})$, $i,j\in{x,y,z}$, some of them necessarily generically vanish or are not independent due to the symmetries of our Hamiltonian (Eq.~3):
\begin{figure}
\includegraphics[width=.6\linewidth]{Fig2.pdf}
\caption{Contribution to the semi-classical response (Eq.~1): In the presence of an illustrative applied field along the $k_y$ direction, the Fermi surface is displaced along the $k_y$-axis in $\mathbf{k}$-space in a manner symmetric about the $k_x=0$ and $k_z=0$ planes. This accounts for the vanishing response $J_x$ and $J_z$.}
    \label{vanishing}
\end{figure}
\begin{equation}
    J_x(A\hat{y})=0,\quad J_z(A\hat{y})=0,\tag{5a}
    \end{equation}

    \begin{equation}
    J_x(A~ \hat{x})=J_z(A\hat{z}),\tag{5b}
    \end{equation}

    \begin{equation}
    J_z(A\hat{x})=J_x(A\hat{z}),\tag{5c}
    \end{equation}

    \begin{equation}
    J_y(A\hat{x})=-J_y(A\hat{z}).\tag{5d}
    \end{equation}
\begin{figure*}
\includegraphics{Fig3.pdf}
\caption{(a) In the presence of an applied transverse field $A\hat{x}$, the Fermi surface of the Hopf link is displaced along the transverse $k_x$-axis in $\mathbf{k}$-space, in turn giving rise to a transverse diagonal optical response $J_x(A\hat{x})$. (b) For the Hopf link, the transverse response weakly depends on $d$, for $r_x=0.8$ and $r_y=1.2>d$. (c) Enhancement of transverse response non-linearity for the Hopf link as $r_x$ increases, for $d=0.5<r_y=1.0$. (d) Enhancement of transverse response non-linearity for the Hopf link with $r_y>d=0.5$ and $r_x=1.0$. (e) Displacement of the Fermi regions of the unlinked nodal loops along the transverse $k_x$ axis in $\mathbf{k}$-space, with the singularities from each loop becoming weaker as the loops are further separated. (f) The transverse response shows a stronger dependence on the $d$ when the nodal loops are no longer topologically linked. As the longitudinal loops separation increases, the response currents become smaller. Here, we contrast the responses of the Hopf link (dashed-dot lines, $d=0.25,0.35$), unlinked nodal loops (dashed lines, $d=1.1,1.3$) and the nodal chain (solid lines $d=0.5,0.9$) for $r_y=0.6$, $r_x=0.4$.
For (b-d), the chosen parameter sets are among those that exhibit the greatest non-linearity in the transverse response. The responses are in units of $env_F$ where $n$ and $v_F$ are the number density and Fermi velocity of the NLSM material.}
\label{Jx}
\end{figure*}

\begin{figure*}
\includegraphics{Fig4.pdf}
\caption{(a) In the presence of an applied longitudinal field $A\hat{y}$, the Fermi surface of the Hopf link is displaced along the longitudinal $k_y$-axis in $\mathbf{k}$-space, in turn giving rise to a longitudinal diagonal optical response $J_y(A\hat{y})$. (b) Enhancement of longitudinal response non-linearity by increasing $d$ in the linked regime with $r_y=1.0>d$, $r_x=0.6$. (c) The shape of the longitudinal response for the Hopf link is independent of the transverse radius $r_x$ but increases in magnitude with $r_x$, as plotted for $r_y=1.2>d=0.8$. (d) Enhancement of longitudinal response non-linearity for the Hopf link by increasing $r_y>d=0.6$ at $r_x=1.0$. (e) Displacement of the Fermi regions of the unlinked nodal loops along the longitudinal $k_y$ axis in $\mathbf{k}$-space, with the singularities from each loop becoming weaker as the loops are further separated.  (f) Comparing the longitudinal responses for the Hopf link (dashed-dot lines $d=0.25, 0.35$), unlinked nodal loops (dashed lines $d=1.2, d=1.4$), and the nodal chain (solid lines $d=0.4, d=0.8$) for $r_y=0.6$, $r_x=0.4$. The longitudinal response decreases in magnitude but becomes more non-linear as $d$ increases. This increasing non-linearity is an artifact of the weaker influence of the singularities as the loops are further separated. For (b-d), the chosen parameter sets are among those that exhibit the greatest non-linearity in the transverse response. The responses are in units of $env_F$ where $n$ and $v_F$ are the number density and Fermi velocity of the NLSM material.
}
\label{Jy}
\end{figure*}
Eq.~5a holds since the nodal loops have mirror symmetry about $k_x=0$ and $k_z=0$ planes respectively. As the dispersion velocity $\langle\mathbf{k}|\hat{\mathbf{J}}|\mathbf{k}\rangle$ (Eq.~1) inherits the same symmetry as the Fermi surface - mirror symmetry about the $k_y$ axis,  When the Fermi surface is displaced along the longitudinal $k_y$ axis (Fig.~\ref{vanishing}), the resulting contribution by the Fermi tubes on the left of the $k_y$ axis will exactly cancel that due to the Fermi tubes on the right. Eq.~5b and Eq.~5c follow from the symmetrical roles played by transverse momenta $k_z$ and $k_z$. Eq.~5d is true because of the loops are related to each other by interchanging $k_x\leftrightarrow-k_z$.  In all, there are only 4 unique currents responses in the following response matrix
\begin{equation}
J_i(A\hat{j})=\begin{pmatrix}J_x(A\hat{x})&J_y(A\hat{x})&J_z(A\hat{x})\\0&J_y(A\hat{y})&0\\J_z(A\hat{x})&-J_y(A\hat{x})&J_x(A\hat{x}).\\
\end{pmatrix}
\end{equation}
Since $A$ scales with $E_0$ for our sinusoidal electric field, Eq. 5 is similar to a conductivity tensor. As we shall show, the diagonal responses are much more affected by the topological linkage of the nodal loops, compared to the non-diagonal response, i.e. Hall responses.

\subsection{Transverse Diagonal Responses $J_x(A\hat{x})$, $J_z(A\hat{z})$}
According to Eq. 1, we generically expect a non-linear semiclassical response since its integral is a highly non-linear function. This is especially the case for a Hopf link, where each nodal loop juts out of the plane of the other nodal loop. As such, it acts as a source of dispersion velocity $\langle\mathbf{k}|\mathbf{\hat{J}}|\mathbf{k}\rangle$, which `destructively interferes' with those due to the nodal loop that encloses this nodal tube~\cite{PhysRevB.102.035138}. When the Fermi surface of the Hopf link starts to translate by the impulse from an applied field (Fig.~\ref{Jx}a), the response, computed from Eq.~1, initially increases sharply followed by a much slower non-monotonic change until the applied field is sufficiently large that the Fermi surface begins to leave the region of influence exerted by the singularity.

Beyond that, the contribution of $\langle\mathbf{k}|\mathbf{\hat{J}}|\mathbf{k}\rangle$ continues to add up, giving a subsequent monotonic response. The extent of non-linearity afforded by each singularity is described by Eq.~2, where the gradient is a sum of the second derivatives of the dispersion. The transverse (perpendicular to the direction of the loop separation) diagonal optical response of the Hopf link was found to observe the following general trends (we plot the most non-linear responses in Figs.~\ref{Jx}b-d):
\begin{enumerate}
    \item The response demonstrates a weak dependence on the longitudinal separation $d$ when the nodal loops are topologically linked,
     but not true otherwise, as illustrated in Figs.~\ref{Jx}b and f respectively.
    \item The extent of response non-linearity is enhanced with larger transverse radius $r_x$ or smaller longitudinal radius $r_y$ i.e. larger loop aspect ratio $r_x/r_y$, as illustrated in Figs.~\ref{Jx}c,d. Correspondingly, the response currents generally decrease with increasing non-linearity.
\end{enumerate}
We can physically understand these trends. Since the relative separation $d$ of the centers of the loops merely changes the relative longitudinal $k_y$ position of the loops, intuitively it should not play a significant role in influencing the transverse response resulting from the transverse translation of the Fermi surface. Indeed, the transverse diagonal response is almost independent of $d$ (Fig.~\ref{Jx}b) in the Hopf link. Yet, this independence from $d$ no longer holds when the loops are no longer topologically linked (Figs.~\ref{Jx}e,f) because further separated loops i.e. sources of singularities correspond to more uniform velocity fields, leading to smaller responses currents. This is unlike the linked cases, where the nodal linkage guarantees the proximity to nodal singularities and topologically ``protects'' the transverse response (previous literature i.e. Ref.~\cite{PhysRevB.102.035138} only reported the topological enhancement of the longitudinal response).

On the other hand, increasing the transverse radius $r_x$ or decreasing the longitudinal radius $r_y$ stretches the aspect ratio of each loop and hence the boundary for ``interference'' (defined by the locus of points where the competition of opposing dispersion velocity vectors terminate) along the transverse $k_x$ direction, giving more room for ``destructive interference'' as the Fermi surface is displaced along the transverse $k_x$ direction. This in turn gives a  greater range of $A\hat{x}$ where the response undergoes a non-monotonic change.

\subsection{Longitudinal Diagonal Response $J_y(A\hat y)$}

Unlike in the transverse diagonal response, the longitudinal separations of the loops $2(r_y-d)$ now greatly influences the longitudinal (along the direction of the loop separation) responses of the Hopf link. Like before, these nodal tubes are sources of `destructive interference' for the dispersion velocity $\langle\mathbf{k}|\mathbf{\hat{J}}|\mathbf{k}\rangle$. The semiclassical response is again computed from Eq. 1 as the Fermi surface of the Hopf link displaces along the longitudinal $k_y$ direction (Fig.~\ref{Jy}a). The general trends for the longitudinal diagonal optical response of the Hopf link are:
\begin{enumerate}
    \item The extent of response non-linearity is enhanced with larger longitudinal loop separation $d$, as illustrated in Fig.~\ref{Jy}b.
		\item Increasing the longitudinal radius $r_y$ also enhances the extent of response non-linearity, as illustrated in Fig.~\ref{Jy}d.
    \item However, the transverse radius $r_x$ does not change the shape of the response curves, but the magnitude of the response does increase with $r_x$, as illustrated in Fig.~\ref{Jy}c.

\end{enumerate}

To physically explain these trends, note that as $d$ increases with $r_y$ fixed, the loop separation $2(r_y-d)$ decreases, enhancing the extent of non-linearity due to the closer proximity of one loop with the band singularity of the other. While the non-linearity is still present, as shown in Fig.~\ref{Jy}f, when the nodal loops are no longer topologically linked, this is an artifact from the weaker dispersion field when the loops are separated farther apart. Increasing $r_y$ at fixed $d$ also enhances non-linear response by providing a greater range of values of $A~\hat{y}$ for significant destructive interference to occur.

On the other hand, $r_x$ only stretches the nodal loop in the transverse $k_x$ direction, which plays no role in the `destructive interference' in the longitudinal direction. But yet, $r_x$ increases the circumferential length of the nodal loop, which gives a greater contribution to $\mathbf{J}$ in Eq.~1 and hence a greater response magnitude.

 \section{Global aspects of optical response anisotropy}
Having discussed the non-linearity of the response tensor along the principal directions, we now present how the anisotropic response behaves as a whole. Due to non-linearity i.e. $\bold J(\bold A_1+\bold A_2)\neq \bold J(\bold A_1)+\bold J(\bold A_2)$, the responses in directions away from the previously studied principal directions $\bold{\hat{A}}=(1,0,0)$, $(0,1,0)$ and $(0,0,1)$ may behave unexpectedly. In the following, we shall represent the response across all directions in the form of constant $|\bold A|=A_0$ level surfaces in $\bold J$-space (recall that $\bold A(t)=\int_{-\infty}^{t}\mathbf{E}(t')dt'$).  We parametrize this surface $\mathcal{S}_A$ with spherical angular coordinates as follows:
\begin{equation}
\mathcal{S}_A=\{A_0\boldsymbol{\hat{A}}=A_0(\cos\theta\cos\phi,\cos\theta\sin\phi,\sin\theta)\}\tag{8}
\end{equation}
where $0\leq\theta\leq\pi,~0\leq\phi\leq 2\pi$. The anisotropy in the optical response of our canonical double nodal loop model (Eqs.~3, 4a,b) due to $\boldsymbol{\alpha}=A_0\boldsymbol{\hat{A}}$ is described by the anisotropy of the smooth $\mathcal{S}_J$ surface embedded in the $J$-space. We represent the azimuthal angle $\phi$ with a colormap defined by $\mathcal{S}_A$ in Fig.~\ref{A0}a, which will be significantly distorted by anisotropy when bijectively mapped onto $\mathcal{S}_J$ (Figs.~\ref{A1} and \ref{A2}).

\begin{figure}
    \centering
 \includegraphics{Fig5.pdf}
    \caption{
    (a) To draw the response surfaces, we represent an external field of constant magnitude $A_0$ but arbitrary direction as a sphere $\mathcal{S}_A$ in $A$-space. The sphere is colored by the azimuthal angle $\phi$, and is significantly distorted by anisotropy when bijectively mapped onto the response surfaces $\mathcal{S}_J$ in 3D $J$-space. (b-c) Illustrative examples of the distortion for our model (Eq.~3) with $r_x=0.6,r_y=1.0,\mu=0.1,A_0=1.25$, for $d=0.5$ and $0.7$ respectively ($\mathbf{J}$ in units of $10^{-3}env_F$). (d) Origin of the thicker ``waist'' of (b) vs. (c) in terms of $J_y(A_0 \hat y)$ response. At $A_0=1.25$, the $d=0.5$ case (b) has a response current almost thrice as large.
     (e) This difference in $J_y(A_0 \hat y)$ corresponds to a large difference in the thickness of the concave ``waist'' of the response surface, as more clearly illustrated by comparing cross sections of (b,c) at $\theta=\frac{\pi}{2}$. }
\label{A0}
\end{figure}

Insight into the shapes of of the constant $A_0$ response surfaces, as well as their significance with regards to nodal topology and geometry, can be obtained by computing the response in a few high-symmetry directions of $\boldsymbol{\hat{A}}$, such as $\boldsymbol{\hat{A}}=\pm(0,1,0),\pm(1,0,1),\pm(1,0,-1)$ (for notational simplicity, we shall henceforth drop the normalization factor). From Eq. 5a, we know that for a longitudinal external field $\boldsymbol{\hat{A}}=\pm(0,1,0)$, there is vanishing transverse response, i.e. $J_x=J_z=0$ for such $\boldsymbol{\hat A}$. Hence the corresponding point on the $\mathcal{S}_J$ surface must be along the $J_y$ axis, with $|J_y|$ being the width of the surface along the $J_y$-axis. This is also the longitudinal diagonal current discussed in Sect. III. Also, the response surface point corresponding to $\boldsymbol{\hat{A}}=\pm(1,0,1)$ must lie in the $J_x-J_z$ plane, since $J_y=0$ from Eq. 5d. And by Eq. 5b, there is an equal contribution of $J_x$ and $J_z$ from $\boldsymbol{\hat{A}}=\pm(1,0,1)$, so all corresponding points from this orientation of $\boldsymbol{\hat A}$ must lie on the 45 degree diagonals along the $J_x$-$J_z$ axes, as evident in Figs.~\ref{A1} and \ref{A2}. Incidentally, points corresponding to $\boldsymbol{\hat A}=\pm(1,0,-1)$ also lie in the $J_x$-$J_z$ plane despite not being directly protected by Eqs.~5a and 5d, a fine illustration of the non-linearity of $\bold J(\bold A)$. Other high symmetry directions which we shall use in the following discussion include $\boldsymbol{\hat{A}}=~\pm(1,1,0),~\pm(0,1,1),~\pm(1,-1,0),~\pm(0,1,-1)$.

Before starting the detailed analysis of the response surface, we consider a quick example. In Sect. III, we highlighted that a greater longitudinal separation $d$ between two topologically linked nodal loops results in a more non-linear response (Fig.~\ref{Jy}b). This is indeed seen in our corresponding $\mathcal{S}_J$ surfaces. At an illustrative $A_0=1.25$, this non-linear response is approximately $3$ times larger i.e. less non-linear for $d=0.5$ vs $d=0.7$ (Fig.~\ref{A0}d), which corresponds to a proportionally smaller width of the ``waist'' of the response surface in the $\hat y$ direction (Figs.~\ref{A0}b,c). This is more evident when viewed in the cross-sectional plane $\hat{A}=(1,0,1)$ with only $\theta=\frac{\pi}{2}$ points plotted (Fig.~\ref{A0}e). That the general shapes of the $d=0.5$ and $d=0.7$ surfaces appear rather similar (Figs.~\ref{A0}b,c) can be understood from the weak dependence on $d$ for the transverse diagonal responses, as discussed in Sect. III. Clearly, the results from Sect. III alone are inadequate in working out the entire response surface and thus the response surfaces provide a bigger picture and as we will see, shed light on nodal topology and geometry.

\begin{figure*}
\includegraphics{Fig6.pdf}
 \caption{
 For the three possible nodal configurations (linked, nodal chain, unlinked) of our double loop Hamiltonian, shown here are their corresponding surfaces $\mathcal{S}_J$ for constant chemical potential $\mu=0.1$ under small, medium and large impulses $A_0=0.1,0.75,1.75$ for the linked and chain case, while $A_0=0.1,0.5,1.0$ for the unlinked case. The qualitative evolution of the response surfaces with increasing $A_0$ differ for the different topological nodal configurations, with the large-$A_0$ response surface elongated towards the $J_z\pm J_x$ direction in the topologically linked/unlinked cases.
 (a-d) The response surfaces for the Hopf link $r_y=0.6>d=0.15$, $r_x=0.6$ exhibit concavity along directions $\pm(0,1,1)$, $(1,0,\pm1)$ due to the locally depressed responses in these directions. (e-h) The response surfaces for the nodal chain $r_y=0.6=d$, $r_x=0.6$ do not obey the same symmetries as those in the Hopf link since the Fermi surface of the nodal chain breaks reflection symmetry. In particular, for intermediate $A_0$ values, the surface looks like a saddle with asymmetric responses, e.g. $\hat{A}=(0,1,1)$ has a greater response than $\hat{A}=(0,-1,-1)$ at $A_0=0.75$. (i-l) The response surfaces for unlinked nodal rings $r_y=0.6>d=0.15$, $r_x=0.6$ do respect these symmetries too but possess topologies distinct to the topologically linked case. The units of the response currents $\mathbf{J}$ are $10^{-3}env_F$, and parameters are chosen such that no periodic images are present. }

\label{A1}
\end{figure*}

\begin{figure*}
\includegraphics{Fig7.pdf}
\caption{The dependence of the surface morphologies on the aspect ratio parameters $r_x,r_y$ of the nodal Hopf link at constant $\mu=0.1$. Again, their response surfaces $\mathcal{S}_J$ are plotted for small, medium and large impulses $A_0=0.1,0.75$ and $1.75$, with current $\mathbf{J}$ are in units of $10^{-3}env_F$. (a-d) The response surfaces for the Hopf link $r_y=0.4>d=0.15$, $r_x=0.6$. (e-h) The same response surfaces as in Fig.\ref{A1}a-d. (i-l) The response surfaces for the Hopf link $r_y=0.6>d=0.15$, $r_x=0.4$ share several similarities with (e-h), but with a more depressed longitudinal response at $A_0=0.75$ and $A_0=1.75$ and a less depressed response along $\hat{A}=\pm(1,0,1)$ at $A_0=1.75$. Overall, the response surfaces differ for different geometric parameters, but do not change as drastically as when the nodal topology changes (Fig.~6).}
\label{A2}
\end{figure*}
\subsection{Response surfaces and nodal topology}
It is already known that~\cite{PhysRevB.102.035138} nodal linkages can enhance the non-linearity of the response, at least in the longitudinal diagonal direction $J_y(A_0\hat{y})$ (Fig.~4). Further, significant non-linearity is also present in the transverse diagonal directions $J_x(A_0\hat{x})$ and $J_z(A_0 \hat{z})$, although not necessarily enhanced by the topological linkage. Hence we shall expect such non-linearity to be manifested in the constant $A_0$ response current surface too.

Indeed, as shown in Fig.~6, typical linked, touching and unlinked nodal loops exhibit significantly different evolutions of the constant $A_0$ response surfaces as $A_0$ is increased. At small $A_0$, all three cases have ellipsoidal-shaped responses surfaces, testimony to anisotropy of the nodal system, even in the linear (small $A_0$) limit. At very large $A_0$ i.e. $A_0=1.75$ where the Fermi regions have been displaced far from their original positions, the response surfaces are all very anisotropic and large, since minimal cancellation of the velocity field $d\varepsilon/d\bold k$ occurs. Their exact shapes depend on the details of the energy dispersion away from the loops, and are non-universal though decidedly anisotropic. What is most interesting is the intermediate $A_0\approx 0.75$ regime, which for the linked case is around when the Fermi region of one loop crosses the singularity from the other loop. For the linked case, the significant non-linearity of the response around intermediate values of $A$ (see Figs.~3 and 4) suppresses the response current, particularly in the longitudinal directions with ``untwisted colors'' in Fig.~7. As such, this leads to a somewhat ``squeezed'' appearance of the response surface compared to that of the nodal chain or unlinked cases, where the response surface looks comparatively ``puffed up'' during the $A_0$ evolution. Compared to the other cases, the unlinked case shows the least variation in response surface shape during the evolution due to the least amount of cancellation of $d\varepsilon/d\bold k$ during the evolution.

A more detailed characterization of the shape of the response surfaces can be performed by analyzing the high symmetry directions. For instance, for the Hopf linked case, the constant $A_0$ surfaces look similar to ellipsoids (Fig.~\ref{A1}b,~\ref{A2}b,f,j) with the longest axis oriented along the 45 degree diagonal between $J_x>0$ and $J_z>0$ axes, corresponding to the direction $\bold J\propto \pm(1,0,1)$. The three `principal axes' of the response surfaces are thus marked out by the three pair of points $\pm(0,1,0)$, $\pm(1,0,1)$ and $\pm(1,0-1)$, as discussed. This is seen in Figs.~\ref{Appendix: ry06rx06response}m-p where the magnitude for the responses along $\pm(1,0,-1)$ are significantly smaller than along $\pm(1,0,1)$ and $\pm(0,1,0)$. This accounts for its oblate appearance.

As we increase $A_0$ for the linked case, the two faces of the constant $A_0$ surface in Fig.~\ref{A1}c which are characterized by $\pm(1,0,-1)$, exhibit concavity, reminiscent of the shape of a red blood cell. The concavity along the directions $\pm(1,0,-1)$ can be better seen in Fig.~\ref{Appendix: diff}b, where the same response surface is viewed from a lateral direction. This characteristic suppressed response is similar to that in Fig.~\ref{A0}d where it was a consequence from the non-monotonicity of the response. In these cases, the responses along $\pm(1,0,-1)$ increase significantly slower than in other directions (Figs.~\ref{Appendix: ry06rx06response}m-p). As $A_0$ continues to increase, the $\hat A$ directions characterized by $\pm(1,0,1)$ and $\pm(0,1,0)$ also exhibit concavity. Again, this is because the responses along $\pm(1,1,0)$, $\pm(0,1-1)$, $\pm(0,1,1)$, $\pm(1,-1,0)$ (which correspond to the eight corners of the surface in Fig.~\ref{A1}d) grow much quicker than along the principal directions, consistent with the individual response curves in Figs.~\ref{Appendix: ry06rx06response}a-l. This surface (Fig.~\ref{A1}d) also show symmetry consistent with the symmetries in the responses (Figs.~\ref{Appendix: ry06rx06response}a-l). For instance, the constant $A_0$ surface has mirror symmetry about the $J_y=0$ plane. Due to the non-negligible thickness of the nodal loops, certain symmetries of the $\mu=0$ nodal system may be broken. This symmetry breaking is particularly pronounced in the nodal chain case, where the two nodal rings intersect with relatively weak dispersion. Generically, the finitely  thick nodal tube can break reflection symmetry in the direction of impulse, i.e. translating the Fermi surface in the $+A_0\hat{y}$ direction results in a different response compared to doing so in the $-A_0\hat{y}$ direction. For small $A_0$, the surfaces are imperfectly rounded blobs (Fig.~\ref{A1}f) and evolve to a saddle shape (Fig.~\ref{A1}g) as $A_0$ grows. This saddle-like shape is more obvious when viewed in a different orientation illustrated in Fig.~\ref{Appendix: diff}c. It can be characterized by 6 out of the 14 points in total: $(0,1,1)$, $(1,1,0)$, $(0,\pm1,0)$, $(0,1,-1)$ and $(-1,1,0)$. Again, we observe asymmetry in the sense that the response grow more than proportionately in the directions $(0,1,1)$, $(1,1,0)$, $(0,1,-1)$ and $(-1,1,0)$ (the negative pair does not grow as fast). As $A_0$ grow, so do the responses corresponding to the remaining 8 points, giving to a surface like in Fig.~\ref{A1}h, with 6 concave sides.

For the unlinked nodal loops, the surface for small $A_0$ is again imperfectly rounded due to linear anisotropy. For intermediate $A_0$, the surface expands into a lemon-like shape due to the relatively small response non-linearity, and can be characterized by the surface directions $\pm(1,1,0),\pm(0,1,1)$ and $\pm(1,-1,0)$. As $A_0$ grow, the responses in the $\pm(1,0,1)$ directions do not grow as fast, and thus the face characterized by this pair of points exhibit concavity. For large $A_0$, the response surface elongates towards the $J_z-J_x$ direction which differs from the linked case where it elongates towards the $J_z+J_x$ direction. This conclusively relates the response surfaces with the topological linkage of the nodal loops.

\subsection{Response surfaces and nodal geometry}

Since the non-linear response current does not correspond to any topologically quantized value, we expect it to be affected by deformations of the nodal structure too. This should apply to both linked and unlinked cases, even if the linked case is more likely to possess strongly non-linear response. Earlier, we have considered circular nodal loops with $r_x=r_y$. How will the surfaces in Fig.~\ref{A1} change as we vary $r_y$ and $r_x$? Varying $r_y$ and $r_x$ add one more layer of complexity in our analysis of these surfaces. To study that, we shall start from the parameters in Fig.~\ref{A1}(a-c) and change $r_y$ and $r_x$ separately, as presented in Fig.~7.

Generally, the surfaces along the same column of Fig.~\ref{A2} will have relatively similar morphologies since they are only slightly distorted from each other and correspond to the same $A_0$ impulse. In particular, the special symmetries in Figs.~\ref{A1}b-d are no longer satisfied as the relative responses along the various directions now behave quite differently. This occurs because a chemical potential of $\mu=0.1$ result in nodal loops with significantly non uniform and asymmetric thickness (Fig.~\ref{A2}i). This further enhances the anisotropy of the response surfaces, particularly at $A_0=0.75$. But at high fields $A_0=1.75$, the surface morphologies of Figs.~\ref{A2}d,h share similar morphological features - suppressed response at $\pm(1,0,1)$ and $\pm(1,0,-1)$. Again, this anisotropic response is due to the varying rates of growth for the responses along the individual directions. For instance, $\pm(1,0,1)$ does not grow as fast as say $\pm(1,1,0)$ in Fig.~\ref{A2}h, but grow at comparable rates in Fig.~\ref{A2}l.

\section{Conclusion}
In this work, we have systematically studied in detail the anisotropic and non-linear optical response of two nodal loops that are linked, unlinked or touching (nodal chain). This system, as parametrized by our canonical two nodal loop model, represents the simplest abstraction of simultaneously occurring nodal loops (linked or unlinked) in nodal materials. First, we studied the effects of nodal geometry and topology individually along various axis directions. Next, we presented constant $A_0$ response surfaces to highlight the anisotropy of the response, and how that global picture can shed light on the overall configuration of the nodal structure. Our findings generalizes existing results on the enhancement of optical response non-linearity by nodal linkages~\cite{PhysRevB.102.035138} to various transverse, Hall and diagonal sectors, and introduces a geometric picture of response non-linearity and anisotropy that will be invaluable in analyzing generic nodal material responses, as well as the engineering of higher-harmonic generation materials for applications like Terahertz radiation generation.

\section{Acknowledgements}
TT is supported by the NSS scholarship by the Agency for Science, Technology and Research (A*STAR), Singapore.
This work is supported by the Singapore Ministry of Education (MOE) Tier 1 start-up grant (R-144-000-435-133).


\section{Introduction}
In the late 19th century, the theory of classical mechanics experienced several issues in reporting the physical phenomena of light masses and high velocity microscopic particles. In 1920s, Bohr’s atomic theory \cite{bohr1928quantum}, Heisenberg’s discovery of quantum mechanics \cite{robertson1929uncertainty} and Schrödinger’s \cite{wessels1979schrodinger} discovery of wave mechanics influence the conception of a new field i.e. the quantum mechanics. In 1982, Feynman \cite{feynman1982simulating} stated that quantum mechanical systems can be simulated by quantum computers in exponential time, i.e. better than with classical
computers. Till then, the concept of quantum computing was thought to be only a theoretical possibility, but over the last three decades the research has evolved such as to make quantum computing applications a realistic possibility \cite{wang2012handbook}.

In the last two decades, the field of swarm intelligence has got overwhelming response
among research communities. It is inspired by nature and aims to build decentralized and self-organized systems by collective behavior of individual agents with each other and with their environment. The research foundation of swarm intelligence is constructed mostly upon two families of optimization algorithms i.e. ant colony optimization   (Dorigo et at. \cite{dorigo1999ant} and Colorni et al. \cite{colorni1992distributed} 1992); and particle swarm optimization (PSO) (Kennedy \& Eberhart \cite{kennedy1995particle} 1995). Originally, the swarm intelligence is inspired by certain natural behaviors of flocks of birds and swarms of ants.

In the mid 1990s, particle swarm optimization technique  was introduced for continuous optimization, motivated by flocking of birds. The evolution of PSO based bio-inspired techniques has been in an expedite development in the last two decades. It has got
attention from different fields such as inventory planning \cite{wang2014modified}, power systems \cite{alrashidi2009survey}, manufacturing       \cite{yildiz2009novel}, communication networks \cite{latiff2007energy}, support vector machines \cite{lin2008particle}, to estimate binary inspiral signal \cite{wang2010particle}, gravitational waves \cite{normandin2018particle} and many more. Similar to evolutionary genetic algorithm, it is inspired by simulation of social behavior, where each individual is called particle, and  group of individuals is called swarm. In multi-dimensional search space, the position and velocity of each particle represent a probable solution. Particles fly around in a search space seeking potential solution. At each iteration, each particle adjusts its position according to the goal of its own and its neighbors. Each particle in a neighborhood share the information with others \cite{sun2004particle}. Later, each particle keeps the record of best solution experienced so far to update their positions and adjust their velocities accordingly.

	\begin{figure*}[!ht]
	\centering
	\includegraphics[scale=1.25]{psoqpso.eps}
	\caption{Particles movement in PSO and QPSO algorithm}
\end{figure*}

Since the first PSO algorithm proposed, the several PSO algorithms have been introduced with plethora of alterations.  Recently, the combination of quantum computing, mathematics and computer science have inspired the creation of optimization techniques. Initially,  Narayanan and Moore \cite{narayanan1996quantum} introduced quantum-inspired genetic algorithm (QGA) in 1995. Later, Sun et al. \cite{sun2004particle} applied the quantum laws of mechanics to PSO and proposed quantum-inspired particle swarm optimization (QPSO). It is the commencement of quantum-behaved optimization algorithms, which has subsequently made a significant impact on the academic and research communities alike.

Recently, Yuanyuan and Xiyu \cite{yuanyuan2018quantum} proposed a quantum evolutionary algorithm to discover communities in complex social networks. Its applicability is tested on five real social networks and results are compared with classical algorithms.  It has been proved that PSO lacks convergence on local optima i.e. it is tough for PSO to come out of the local optimum once it confines into optimal local region. QPSO with mutation operator (QPSO-MO) is proposed to enhance the diversity to escape from local optimum in search \cite{liu2005quantum}. Protopopescu  and Barhen \cite{protopopescu2002solving} solved set of global optimization problems efficiently using quantum algorithms. In future, the proposed algorithm can be integrated with matrix product state based quantum classifier  for supervised learning \cite{40, 44}.

In this paper, we have combined QPSO with Cauchy mutation operator to add long jump ability for global search and natural selection mechanism for elimination of particles. The results shown that it has great tendency to overcome the problem of trapping into local search space. Therefore, the proposed hybrid QPSO strengthened the local and global search ability and outperformed the other variants of QPSO and PSO due to fast convergence feature.

The illustration of particles movement in PSO and QPSO algorithm is shown in Fig 1. The big circle at center denotes the particle with the global position and other circles are particles. The particles located away from global position are lagged particles.  The blue color arrows signify the directions of other particles and the big red arrows point towards the side in which it goes with high probability. During iterations, if the lagged particle is unable to find better position as compared to present global position in PSO, then their impact is null on the other particles. But, in QPSO, the lagged particles move with higher probability in the direction of gbest position. Thus, the contribution of lagged particles is more to the solution in QPSO in comparison with PSO algorithm.

The organization of rest of this paper is as follows: Sect. 2 is devoted to prior work. In Sect. 3, the quantum particle swarm optimization is described. In
Sect. 4, the proposed hybrid QPSO algorithm with Cauchy distribution and natural selection mechanism  is presented. The experimental results are
plotted for a set of benchmark problems and compared with several QPSO variants in Sect. 5. The correctness and time complexity are analyzed in Section 6. QPSO-CD is applied to three constrained engineering design problems in Sect. 7.
Finally, Sect. 8 is the conclusion.

\section{Prior Work}
Since the quantum-behaved particle swarm optimization was proposed, various revised variants have been emerged. Initially, Sun et al. \cite{sun2004particle} applied the concept of quantum computing to PSO and developed a quantum Delta potential well model for classical PSO \cite{sun2004global}. It has been shown that the convergence and performance of QPSO are superior as compared to classical PSO. The selection and control of parameters can improve its performance, which is posed as an open problem. Sun et al. \cite{sun2007using} tested the performance of QPSO on constrained and unconstrained problems. It has been claimed that QPSO is a promising optimization algorithm, which performs better than classical PSO algorithms. In 2011, Sun et al. \cite{sun2011quantum} proposed QPSO with Gaussian distribution (GAQPSO) with the local attenuator point and compared its results with several PSO and QPSO counterparts. It has been proved that GAQPSO is efficient and stable with superior features in quality and robustness of solutions.

Further, Coelho \cite{dos2010gaussian} applied GQPSO to constrained engineering problems and shown that the simulation results of GQPSO are much closer to the perfect solution with small standard deviation. Li et al. \cite{li2012improved} presented a cooperative QPSO using Monte Carlo method (CQPSO), where particles cooperate with each other to enhance the performance of original algorithm. It is implemented on several representative functions and performed better than the other QPSO algorithms in context of computational cost and quality of solutions. Peng et al. introduced \cite{peng2013quantum} QPSO with Levy probability distribution and claimed that there are very less chances to  be stuck in local optimum.

Researchers have applied PSO and QPSO to real-life problems and achieved optimal solutions as compared to existing algorithms. Ali  et al. \cite{90} performed energy-efficient clustering in mobile ad-hoc networks (MANET) with PSO. The similar approach can be followed to analyse and execute mobility over MANET with QPSO-CD \cite{88}. Zhisheng \cite{zhisheng2010quantum} used QPSO in economic load dispatch for power system and proved superior to other existing PSO optimization algorithms. Sun et al. \cite{sun2006qpso} applied QPSO for QoS multicast routing. Firstly, the QoS multicast routing is converted into constrained integer problems and then effectively solved by QPSO with loop deletion task. Further, the performance is investigated on random network topologies. It has been proved that QPSO is more powerful than PSO and genetic algorithm.
Geis and Middendorf proposed PSO with Helix structure for finding ribonucleic acid (RNA) secondary structures with same structure and low energy \cite{geis2011particle}. The QPSO-CD algorithm can be used with two-way quantum finite automata to model the RNA secondary structures \cite{bhatia2018modeling}.
 Bagheri et al. \cite{bagheri2014financial} applied the QPSO for tuning the parameters of adaptive network-based fuzzy inference system (ANFIS) for forecasting the financial prices of future market. Davoodi et al. \cite{davoodi2014hybrid} introduced a hybrid improved QPSO with Neldar Mead simplex method (IQPSO-NM), where NM method is used for tuning purpose of solutions. Further, the proposed
algorithm is applied to solve load flow problems of power system and acquired the convergence accurately with efficient search ability. Omkar \cite{omkar2009quantum} proposed QPSO for multi-objective design problems and results are compared with PSO. Recently, Fatemeh et al. \cite{fatemeh2019shuffled} proposed QPSO with shuffled complex evolution (SP-QPSO) and its performance is demonstrated using five engineering design problems. Prithi and Sumathi \cite{prithi2020ld2fa} integrated the concept of classical PSO with deterministic finite automata for transmission of data and intrusion detection. The proposed algorithm QPSO-CD can be used with quantum computational models for wireless communication \cite{10, 20, 30}.

\section{Quantum Particle Swarm Optimization}
Before we explain our hybrid QPSO-CD algorithm mutated with Cauchy operator and natural selection method, it is useful to define the notion of quantum PSO. We assume that the reader is familiar with the concept of classical PSO; otherwise, reader can refer to particle swarm optimization algorithm \cite{kennedy2010particle, shi2001particle}. The specific principle of quantum PSO is given as:

In QPSO, the state of a particle can be represented using wave function $\psi(x, t)$. The probability density function $|\psi(x, t)|^2$ is used to determine the probability of particle occurring in  position \textit{x} at any time \textit{t} \cite{sun2004particle, sun2006qpso}. The position of particles is  updated according to equations:

\begin{equation}
x_{i, j}(t+1)=p_{i, j} (t) \pm \alpha . |mbest_{i, j}(t)-x_{i, j}(t)|. ln(1/u)
\end{equation}
\begin{equation}
p_{i, j}(t)=  (\phi.P_{i, j}(t)+(1-\phi).G_j(t)),  (1 \leq i \leq  N, 1 \leq j \leq M)
\end{equation}
where each particle must converge to its local attractor $p=(p_1, p_2,..., p_D)$, where \textit{D} is the dimension, \textit{N} and \textit{M} are the number of particles and iterations respectively, $P_{i, j}$ and $G_{j}$ denote the previous and optimal position vector of each particle respectively, $\phi=c_1.r_1/(c_1r_1+c_2r_2)$, where $c_1$; $c_2$ are the acceleration coefficients, $r_1$; $r_2$ and \textit{u} are normally distributed random numbers in (0, 1), $\alpha$ is contraction-expansion coefficient and mbest defines the mean of best positions of particles as:
\begin{equation}
mbest_{i, j}(t)=\dfrac{1}{N} \sum_{i=1, j=1}^{N, M}P_{i, j}(t)=\bigg(\dfrac{1}{N}\sum_{i=1}^{N}P_{i, 1}(t), \dfrac{1}{N}\sum_{i=1}^{N}P_{i, 2}(t), . . ., \dfrac{1}{N}\sum_{i=1}^{N}P_{i, D}(t)  \bigg)
\end{equation}
In Eq. (1), $\alpha$ denotes contraction-expansion coefficient, which is setup manually to control the speed of convergence. It can be decreased linearly or fixed. In PSO, $\alpha < 1.782$ to ensure convergence performance of the particle. In QPSO-CD, the value of $\alpha$ is determined  by $\alpha$=1-(1.0-0.5)\textit{k}/\textit{M}, i.e. decreases linearly from 1.0 to 0.5 to attain good performance, where \textit{k} is present iteration and \textit{M} is maximum number of iterations.

\section{Hybrid Particle Swarm Optimization}
The hybrid quantum-behaved PSO algorithm with Cauchy distribution and natural selection strategy (QPSO-CD) is described as follows:

The QPSO-CD algorithm begins with the standard QPSO using equations (1), (2) and (3). The position and velocity of particles cannot be determined exactly due to varying dynamic behavior.  So, it can only be learned with the probability density function. Each particle can be mutated with Gaussian or Cauchy distribution. We mutated QPSO with Cauchy operator due to its ability to make larger perturbation.  Therefore, there is a higher probability with Cauchy as compared to Gaussian distribution to come out of the local optima region. The QPSO algorithm is mutated with Cauchy distribution to increase its diversity, where mbest or global best position is mutated with fixed mutation probability ({\textit{Pr}). The probability density function $(f(x))$ of the standard Cauchy distribution is given as:

	\begin{equation}
	f(x)=\dfrac{1}{\pi(1+x^2)} ~ ~ ~ ~ -\infty < x < \infty,
	\end{equation}

It should be noted that mutation operation is executed on each vector by adding Cauchy distribution random value (\textit{D}(.)) independently such that
\begin{equation}
x^{'}=x+\phi D(.)
\end{equation}
where $x^{'}$ is new location after mutated with random value to \textit{x}. At last, the position of particle is selected and the particles of swarm are sorted on the basis of their fitness values after each iteration. Further, substitute the group of particles having worst fitness values with the best ones and optimal solution is determined.  The main objective of using natural mechanism is to refine the capability and accuracy of QPSO algorithm.

 \begin{center}
\begin{minipage}{.7\linewidth}
\begin{algorithm}[H]
		\caption{QPSO-CD algorithm}
		\label{pseudoPSO}
		\begin{algorithmic}[1]
			\State The swarm is initialized with random numbers distributed uniformly: random $x_i$.
			\State Do
			\State $\alpha$ decreases linearly from 1.0 to 0.5

			\For {\textit{k}=1 to \textit{M} do}
			\State $\alpha=1-(1.0-0.5).k/M$
			\If {$Pr<rand(0,1)$}
			\State Calculate the mbest of the swarm using Eq (3)
			\EndIf
			\For {\textit{i}=1: to \textit{N} do}
			\If {Fitness$(p_i) <$ Fitness$(x_i)$} $p_i=x_i$;
			\State $G=argmin(Fitness(p_i)$;
			\EndIf
			\For{\textit{j}:1 to \textit{D} do}
			\State $r_{1}=rand(0, 1); r_{2}=rand(0,1)$;
			\State $\phi=c_{1}.r_{1}/(c_{1}r_{1}+c_{2}r_{2})$;
			\State $p_{i, j}=\phi.P_{i, j}+(1-\phi). G_j$;
			\If {$rand(0,1) < 0.5$}
			\State $x_{i,j}=p_{i, j}+\alpha. abs(mbest_{i, j}-x_{i, j}). log (1/u)$
			\Else
			\State $x_{i,j}=p_{i, j}-\alpha. abs(mbest_{i, j}-x_{i, j}). log (1/u)$
			\EndIf
			\EndFor
			\State $Fx(i)$=Fitness$(x(i,:))$;
			\EndFor

			\State $[SF, Sx]=sort(Fx)$;
			\State $Z=round((N-1)/S)$; \Comment{\textit{S} is selection parameter}
			\State $x(Sx((N-Z+1):N))=x(Sx(1:Z))$;  \Comment{Sort the particles from best to worst position}
			\EndFor
			\State Until termination criterion is met

		\end{algorithmic}
	\end{algorithm}
	 \end{minipage}
	 \end{center}

The natural selection method is used to enhance the convergence characteristics of proposed QPSO-CD algorithm, where the fitter solutions are used for the next iteration. The procedure of selection method for \textit{N} particles is as follows:
\begin{equation}
F(X(t))=\{F(x_1(t)), F(x_2(t)), ..., F(x_N(t))\}
\end{equation}
where $X(t)$ is position vector of particles at time \textit{t} and $F(X(t))$ is the fitness function of swarm.  Next step is to sort the particles according to their fitness values from best one to worst position such that
\begin{equation}
\begin{split}
F(X^{'}(t))=\{F(x_1^{'}(t)),   F(x_2{'}(t)) , ..., &F(x_N{'}(t))\} \\
X^{'}(t)=\{x_1^{'}(t), x_2{'}(t), ..., x_N{'}(t)\}
\end{split}
\end{equation}

In Algorithm 1, $SF$ and $Sx$ are the sorting functions of fitness and position respectively. On the basis of natural selection parameters and fitness values, the positions of swarm particles are updated for the next iteration,
\begin{equation}
\begin{split}
X^{'}(t)=\{x_1^{''}(t), x_2{''}(t), ...,  x_S{''}(t)\}, & \\
X^{''}_{k}(t)=\{x_1^{'}(t), x_2{'}(t), ..., x_Z{'}(t)\}
\end{split}
\end{equation}
where $(1 \leq k \leq S)$, \textit{S} denotes the selection parameter, \textit{Z} signifies the number of best positions selected according to fitness values such that $S=N/Z$ and $X^{''}(t)$ is updated position vector of particles. The selection parameter \textit{S} is generally set as 2 to replace the half of worst positions with the half of best positions of particles. It improves the precision of the direction of particles, protect the global searching capability and speed up the convergence.

\section{Experimental results}

The performance of proposed QPSO-CD algorithm is investigated on representative benchmark functions, given in Table 1. Further, the results are compared with classical PSO (PSO), standard QPSO, QPSO with delta potential (QDPSO) and QPSO with mutation operator (QPSO-MO). The details of numerical  benchmark functions are given in Table 1.
\begin{table}[!ht]
	\centering
	\caption{Details of benchmark functions}
	\begin{tabular}{|c|c|}
		\hline
		Test function  & Initial range  \\
		\hline
		\multicolumn{2}{|c|}{\textbf{Sphere function}}\\
		\hline
		$f_1(x)=\sum\limits_{i=1}^{n}x_{i}^2$ & (-100, 100)  \\
		\hline
		\multicolumn{2}{|c|}{\textbf{Rosenbrock function}} \\
		\hline
		$f_2(x)=\sum\limits_{i=1}^{n}100(x_{i+1}-x_{i}^2)^2+(x_{i}-1)^2$ & (-5.12, 5.12)  \\
		\hline
		\multicolumn{2}{|c|}{\textbf{Greiwank function}}\\
		\hline
		$f_3(x)=\dfrac{1}{4000} \sum\limits_{i=1}^{n}x_{i}^2-\prod_{i=1}^{n}cos(\dfrac{x_i}{\sqrt{i+1}})+1$& (-600, 600) \\
		\hline
		\multicolumn{2}{|c|}{\textbf{Rastrigrin function}} \\
		\hline
		$f_4(x)=\sum\limits_{i=1}^{n}(x_{i}^2-10cos(2\pi x_{i})+10)$ & (-5.12, 5.12) \\
		\hline
	\end{tabular}
\end{table}

The performance of QPSO has been widely tested for various test functions. Initially, we have considered four representative benchmark functions to determine the reliability of QPSO-CD algorithm. For all the experiments, the size of population is 20, 40 and 80 and dimension sizes are 10, 20 and 30. The parameters for QPSO-CD algorithm are as follows: the value of $\alpha$ decreases from 1.0 to 0.5 linearly, the natural selection parameter \textit{S}=2 is taken, $c_1, c_2$ correlation coefficients are set equal to 2.

The mean best fitness values of PSO, QPSO, QDPSO, QPSO-MO and QPSO-CD are recorded for 1000, 1500 and 2000 runs of each function. Fig. 2 to Fig. 5 depict the performance of functions $f_1$ to $f_4$ with respect to mean best fitness against the number of iterations. In Table 2, \textit{P} denotes the population, dimension is represented by \textit{D} and \textit{G} stands for generation. The numerical results of QPSO-CD shown optimal solution with fast convergence speed and high accuracy. The results shown that QPSO-CD performs better on Rosenbrock function than its counterparts in some cases. When the size of population is 20 and dimension is 30, the results of proposed algorithm are not better than QPSO-MO, but QPSO-CD performs better than PSO, QPSO and QDPSO. The performance of QPSO-CD is significantly better than its variants on Greiwank and Rastrigrin functions. It has outperformed other algorithms and obtained optimal solution (near zero) for Greiwank function. In most of the cases, QPSO-CD is more efficient and outperformed the other algorithms.

\begin{figure}
\centering
\begin{minipage}{.5\textwidth}
  \centering
  \includegraphics[width=1\linewidth]{f1.eps}
  \captionof{figure}{Effectiveness of QPSO-CD for sphere function $f_1$}
  \label{fig:test1}
\end{minipage}%
\begin{minipage}{.5\textwidth}
  \centering
  \includegraphics[width=1\linewidth]{f2.eps}
  \captionof{figure}{Effectiveness of QPSO-CD for Rosenbrock function $f_2$}
  \label{fig:test2}
\end{minipage}
\end{figure}

\begin{figure} [!ht]
\centering
\begin{minipage}{.5\textwidth}
  \centering
  \includegraphics[width=1\linewidth]{f3.eps}
  \captionof{figure}{Effectiveness of QPSO-CD for Greiwank function $f_3$}
  \label{fig:test1}
\end{minipage}%
\begin{minipage}{.5\textwidth}
  \centering
  \includegraphics[width=1\linewidth]{f4.eps}
  \captionof{figure}{Effectiveness of QPSO-CD for Rastrigrin function $f_4$}
  \label{fig:test2}
\end{minipage}
\end{figure}
\begin{table}[!ht]
	\centering
	\caption{Comparison results of Sphere and Rosenbrock functions}
	\begin{tabular}{c|c|c|c|c|c|c|c|c|c|c|c|c}
		\hline
		\multicolumn{3}{c|}{}& \multicolumn{5}{c|}{\textbf{Sphere function}} &  \multicolumn{5}{c}{\textbf{Rosenbrock function}}\\
		\hline
		P & D & G & PSO & QPSO & QDPSO & QPSO-MO & QPSO-CD & PSO & QPSO & QDPSO & QPSO-MO & QPSO-CD \\
		\hline
		\multirow{3}{*}{20}& 10 & 1000 &0.0  & 4.01e-40  & 1.513e-49  &1.508e-48  & 1.738e-50 & 95.10 & 58.41 & 14.22 & 22.18 & 34.67  \\
		& 20 & 1500 & 0.0 & 2.58e-21  & 1.339e-30   & 1.296e-31 & 1.032e-30 & 204.38 & 110.5 & 175.31 & 68.40 & 54.76 \\
		& 30 & 2000 & 0.0 & 2.08e-13  & 1.953e-21  & 1.918e-21  & 1.808e-21 & 314.46 & 148.5 & 242.37 & 113.30 & 122.5 \\
		\hline
		\multirow{3}{*}{40}& 10 & 1000 & 0.0 & 2.73e-67 & 1.087e-73  & 1.146e-51 & 1.154e-72 & 70.28 & 10.42 & 15.86 & 7.985 &  8.843    \\
		& 20 & 1500 & 0.0  & 4.84e-28 & 1.397e-42  & 1.417e-42  & 1.237e-41 & 178.98 & 48.45 & 112.46 & 52.93 & 41.77 \\
		& 30 & 2000 & 0.0 & 2.02e-25  & 2.850e-30  &  2.471e-28 & 1.946e-23  & 288.58 & 58.32 & 76.42 & 64.19& 58.04 \\
		\hline
		\multirow{3}{*}{80}& 10 & 1000 &  0.0 & 7.66e-95 & 5.553e-90  & 4.872e-71 &  6.437e-72 & 36.29 & 8.853 & 36.34 & 5.715 & 7.419  \\
		& 20 & 1500 & 0.0 & 1.62e-60 &  1.654e-54  & 1.677e-58  & 1.609e-62 & 84.78 & 34.88 & 23.54 &24.45 & 21.78 \\
		& 30 & 2000 & 0.0 & 2.05e-44 & 1.042e-40  & 1.131e-42  & 1.128e-41 & 202.58 & 52.17 & 70.81&45.22 & 40.97 \\
		\hline
	\end{tabular}
\end{table}

\begin{table}[!ht]
	\centering
	\caption{Comparison results of Greiwank and Rastrigrin functions}
	\begin{tabular}{c|c|c|c|c|c|c|c|c|c|c|c|c}
		\hline
		\multicolumn{3}{c|}{}& \multicolumn{5}{c|}{\textbf{Greiwank function}} &  \multicolumn{5}{c}{\textbf{Rastrigrin function}}\\
		\hline
		P & D & G & PSO & QPSO & QDPSO & QPSO-MO & QPSO-CD & PSO & QPSO & QDPSO & QPSO-MO & QPSO-CD \\
		\hline
		\multirow{3}{*}{20}& 10 & 1000 & 0.089 & 0.078 & 0.1003 &0.0732 & 0.072 & 5.526 & 5.349 & 4.969 & 4.478 & 4.051   \\
		& 20 & 1500 & 0.0300 & 0.2001 &  0.0086 & 0.0189 & 0.0078 & 23.17 & 21.28 & 17.08 & 15.63 & 13.22\\
		& 30 & 2000 & 0.0181 & 0.0122 & 0.0544& 0.0103 & 0.0026 & 46.29 & 32.57 & 48.61 & 27.80 & 31.48 \\
		\hline
		\multirow{3}{*}{40}& 10 & 1000 &0.0826 & 0.055 & 0.048 & 0.0520 &  0.041 & 3.865 & 3.673 & 2.032 & 3.383 &  2.100 \\
		& 20 & 1500 & 0.0272 & 0.0149 & 0.0004 & 0.0247& 0.0106 & 15.68 & 14.37 & 10.94 & 11.01 & 10.77\\
		& 30 & 2000 & 0.0125 & 0.0117 & 0.0009& 0.0105& 0.0102 & 37.13 & 23.01 & 21.37& 21.01& 21.19\\
		\hline
		\multirow{3}{*}{80}& 10 & 1000 & 0.0723 & 0.0341 &0.0 & 0.0542& 0.0702 & 2.562 & 2.234 & 0.923 & 2.183&  1.943  \\
		& 20 & 1500 & 0.0274  &  0.0189 & 0.0 & 0.0194& 0.0161 & 12.35 & 9.66 & 6.955 & 8.075& 7.021 \\
		& 30 & 2000 & 0.0123 & 0.0118 & 0.0 & 0.0082 & 0.0031 & 26.89 & 17.48 & 18.13 &14.99 & 11.73 \\

		\hline
	\end{tabular}
\end{table}

\section{Correctness and Time Complexity Analysis of a QPSO-CD Algorithm}
In this Section, the correctness and time complexity of a proposed algorithm QPSO-CD is analyzed and compared with the classical PSO algorithm.
\begin{theorem}
	The sequence of random variables $\{S_n, n \geq 0\}$ generated by QPSO with Cauchy distribution converges to zero in probability as \textit{\textit{n}} approaches infinity.
		\end{theorem}
	\begin{proof}
		Recall, the probability density function of standard Cauchy distribution and its convergence probability \cite{rudolph1997local} are given as
	\begin{equation}
	\begin{gathered}
	f(s)=\dfrac{1}{\pi (1+s^2)} ~ \text{for} -\infty < s < \infty, \\
P(x \leq S_n \leq y) = \dfrac{1}{\pi} \int_{y}^{x} \dfrac{ds}{(1+s^2)}, ~ \forall x \leq y
\end{gathered}
	\end{equation}
		Consider a random variable $Q_n$ interpreted as
		$$ Q_n = \alpha S_n, ~ \alpha =\dfrac{1}{n^\lambda} $$
		where $\lambda$ denotes a fixed positive constant. Correspondingly, the probability density function can be calculated as
		\begin{align*}
		P(Q_n \leq q) & = P(\alpha S_n \leq q) \\
		& = P\left(S_n \leq \dfrac{q}{\alpha}\right)\\
		& = \int_{-\infty}^{\dfrac{q}{\alpha}} f(s) ds. \dfrac{dP(q_n \leq q)}{dq}\\
		& = \dfrac{1}{\alpha} f\left(\dfrac{q}{\alpha}\right)
		\end{align*}
	i.e. the probability density function of random variable $Q_n$.
	\begin{align*}
	P(|Q_n|>\xi) & = P(|S_n|>\xi n^\lambda)\\
	& = P(S_n > \xi n^\lambda)+ P(S_n > - \xi n ^\lambda)\\
	& = P(\xi n^\lambda < S_n < \infty)+ P(- \infty < S_n < - \xi n^\lambda)
   \end{align*}
		Using Eq. (9), the probability density function of random variable $Q_n$ becomes
		\begin{align*}
			P(|Q_n|>\xi) & = \dfrac{1}{\pi}\int_{\xi n^\lambda}^{\infty}\dfrac{ds}{\pi (1+s^2)} + \dfrac{1}{\pi}\int_{\infty}^{-\xi n^\lambda}\dfrac{ds}{\pi (1+s^2)} \\
			& = \left[1+\dfrac{1}{\pi} \int_{\xi n^\lambda}^{-\xi n^\lambda}\dfrac{ds}{\pi (1+s^2)}\right]\\
			& = 1-\dfrac{1}{\pi} \int_{-\xi n^\lambda}^{\xi n^\lambda}\dfrac{ds}{\pi (1+s^2)}= 0 ~ \text{as} ~ n \rightarrow \infty
		\end{align*}

		This completes the proof of the theorem.
	\end{proof}
\begin{defn}
	Let $\{S_n\}$ a random sequence of variables. It is converges to some random variable \textit{s} with probability 1, if for every $\xi > 0$ and $\lambda >0$, there exists $n_1(\xi, \lambda)$ such that $P(|S_n-s|< \xi) > 1-\lambda, \forall n > n_1$ or
	\begin{equation}
	P\left(\lim_{n \rightarrow \infty} |S_n-s|< \xi\right)=1
	\end{equation}
\end{defn}

The efficiency of the QPSO-CD algorithm is evaluated by number of steps needed to reach the optimal region $R(\xi)$. The method is to evaluate the distribution of number of steps needed to hit $R(\xi)$ by comparing the expected value and moments of distribution. The total number of stages to reach the optimal region is determined as $W(\xi)=inf\{n \mid f_n  \in R(\xi)\}$. The variance $V(W(\xi))$ and expectation value $E(W(\xi))$ are determined as

\begin{equation}
E(W(\xi))= \sum_{n=0}^{\infty}nx_n,
\end{equation}

\begin{equation}
\begin{split}
V(W(\xi)) & =E(W^2(\xi))-\{E(W(\xi))\}^2\\
& = \sum_{n=0}^{\infty}n^2x_n-\left(\sum_{n=0}^{\infty}nx_n\right)^2
\end{split}
\end{equation}

In fact, the $E(W(\xi))$ depends upon the convergence of $\sum_{n=0}^{\infty}nx_n$. It is needed that $\sum_{j=0}^{\infty}x_n=1$, so that QPSO-CD can converge globally. The number of objective function  evaluations are used to measure time. The main benefit of this approach is that it shows relationship  between processor and measure time as the complexity of objective function increases. We used Sphere function  $f(x)=x^{T}.x$ with a linear constraint $g(x)=\sum^{n}_{j=0}x_j \geq 0$ to compute the time complexity. It has minimum value at 0. The value of optimal region is set as $R(\xi)=R(10^{-4}).$ To determine the time complexity, the algorithms PSO and QPSO-CD are executed 40 times on $f(x)$ with initial scope [-10, 10]$^{N}$, where \textit{N} denotes the dimension. We determine the mean number of objective function evaluations ($W(\xi)$), the variance ($V(W(\xi))$), the standard deviation (SD) ($\sigma_{W(\xi)}$), the standard error (SE) ($\sigma_{W(\xi)}/ \sqrt{40}$)  and ratio of mean and dimension ($W(\xi)/N$). The contraction coefficient $\alpha= 0.75$ is used for QPSO-CD and constriction coefficient $\chi=0.73$ for PSO with acceleration factors $c_1=c_2$=2.25.

\begin{table}[!ht]
	\centering
	\caption{Results of the time complexity for QPSO-CD algorithm}
	\begin{tabular}{c|c|c|c|c|c}
		\hline
		Dimension (N) & Mean & Variance & SD & SE & Mean/N  \\
		\hline
	2	 &302.38  & 4164.3 & 64.532   &10.203 & 151.19 \\
	3	& 452.18   & 4541.9 & 67.394   &10.655 & 150.72\\
	 4 & 621.29   & 5208.2 & 72.168 &  11.410  & 155.32\\
	 5 &755.88 & 6675.0 & 81.701   & 12.918 & 151.17\\
	 6 & 879.13 & 8523.7 & 92.324    & 14.597 & 146.52\\
	 7 &  1022.06 & 9575.4  & 97.854      & 15.472 & 146.00\\
	 	8 & 1158.52  & 10269.7  & 101.341   & 16.023 & 144.81\\
	 	9 & 1308.17  & 12053.4 &   109.788  & 17.359 & 145.35 \\
	 	10 & 1459.3 & 12648.3  &  112.465   & 17.782 & 145.93\\
		\hline
	\end{tabular}
\end{table}

\begin{table}[!ht]
	\centering
	\caption{Results of the time complexity for PSO algorithm}
	\begin{tabular}{c|c|c|c|c|c}
		\hline
		Dimension (N) & Mean & Variance & SD & SE & Mean/N  \\
		\hline
		2	 &691.4  & 17297.5 & 131.52   & 20.795 &  345.7\\
		3	& 979.1   & 22281.5&  149.27  & 23.601& 326.3\\
		4 & 1167.2   & 24282.9 & 155.72 &  24.638 & 291.8\\
		5 & 1328.7 & 21853.7 & 147.83  &  23.373& 265.7\\
		6 & 1489.9 &  32008.7&   178.91 &  28.288& 248.3 \\
		7 &  1744.3 & 502297  &  224.12    & 35.436 & 249.1 \\
		8 & 1978.5  &  41233.3 &  203.06  & 31.106 & 247.3 \\
		9 & 2259.1  & 36217.8 & 190.31  & 30.090&  251.0\\
		10 & 2604.2 &  43559.8 &  208.71   &  32.999& 260.4\\
		\hline
	\end{tabular}
\end{table}
\begin{figure*}[!ht]
	\centering
	\includegraphics[scale=0.6]{Time.eps}
	\caption{Time complexity results for PSO and QPSO-CD}
\end{figure*}
\begin{figure*}[!ht]
	\centering
	\includegraphics[scale=0.6]{Roceffi.eps}
	\caption{Comparison of Correlation coefficient of PSP and QPSO-CD}
\end{figure*}

Table 1 and 2 show the statistical results of time complexity test for QPSO-CD and PSO algorithm, respectively. Fig 6 indicates that the time complexity of proposed algorithm increases non-linearly as the dimension increases. However, the time complexity of PSO algorithm increases adequately linearly. Thus, the time complexity of QPSO-CD is lower than PSO algorithm. In Fig 7, QPSO-CD shows a strong correlation between $W(\xi)$ and \textit{N}, i.e. the correlation coefficient is 0.9996. For PSO, the linear correlation coefficient is 0.9939, which is not so phenomenal as that in case of QPSO-CD.  The relationship between  mean and dimension clearly shows that the value of correlation coefficient is fairly stable for QPSO-CD as compared to PSO algorithm.
\section{QPSO-CD for Constraint Engineering Design Problems}
 There exists several approaches for handling constrained optimization problems. The basic principle is to convert the constrained optimization problem to unconstrained by combining objective function and penalty function approach. Further, minimize the newly formed objective function with any unconstrained algorithm. Generally, the constrained optimization problem can be described as in Eq. (13).

The objective is to minimize the objective function \textit{f(x)} subjected to equality $(h_j(x))$ and inequality $(g_i(x))$ constrained functions, where $p(i)$ is the upper bound and $q(i)$ denotes the search space lower bound. The strict inequalities of form $g_i(x) \geq 0$ can be converted into $-g_i(x)\leq 0$ and  $h_i(x)$ equality constraints can be converted into inequality constraints $h_i(x) \geq 0$ and $h_i(x) \leq 0$. Sun et al. \cite{sun2007using} adopted non-stationary penalty function to address non-linear programming problems using QPSO. Coelho \cite{dos2010gaussian} used penalty function with some positive constant i.e. set to 5000. We adopted the same approach and replace the constant with dynamically allocated penalty value.
\begin{equation}
\begin{split}
& \min\limits_{x}=f(x) \\
\text{subject to}  \\
& g_i(x) \leq 0, i=0, 1,... n-1\\
& h_j(x)=0, j=1,2,...r \\
& p(i)\leq x_i \leq q(i), 1 \leq i \leq m\\
& x=\{x_1,x_2, x_3,..., x_m\}
\end{split}
\end{equation}

Usually, the procedure is to find the solution for design variables that lie in search space upper and lower bound constraints such that $x_i \in [p(i), q(i)]$. If solution violates any of the constraint, then the following rules are applied
\begin{equation}
\begin{split}
& x_i=x_i + \{p(x_i)-q(x_i)\}.~ rand[0,1] \\
& x_i=x_i - \{p(x_i)-q(x_i)\}.~ rand[0,1]
\end{split}
\end{equation}

where rand[0, 1] is randomly distributed function to select value between 0 and 1. Finally, the unconstrained optimization problem is solved using dynamically modified penalty values according to inequality constraints $g_i(x)$. Thus, the objective function is evaluated as

\begin{equation}
F(x)=\left\{
\begin{array}{ll}
f(x)~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~\text{if}~ g_i(x)\leq 0\\
f(x)+y(t). \sum\limits_{i=1}^{n} g_i(x)~ ~  \text{if}~ g_i(x)>0\\
\end{array} \right \}
\end{equation}

where $f(x)$ is the main objective function of optimization problem in Eq. (13), \textit{t} is the iteration number and $y(t)$ represents the dynamically allocated penalty value.

In this Section, QPSO-CD is tested for three-bar truss, tension/compression spring  and pressure vessel design problems consisting
different members and constraints. The performance of QPSO-CD is compared and analyzed with the results of PSO, QPSO, and SP-QPSO algorithms as reported in the literature.
\subsection{Three-bar truss design problem}
Three-bar truss is a constraint design optimization problem, which has been widely used to test several methods. It consists cross-section areas of three bars $x_1$ (and $x_3$) and $x_2$ as design variables. The aim of this problem is to  minimize the weight of truss subject to maximize the stress on these bars. The structure should be symmetric and subjected to two constant loadings $P_1=P_2=P$ as shown in Fig 8. The mathematical formulation of two design bars ($x_1$, $x_2$) and three restrictive mathematical functions are described as:

\begin{equation}
\begin{split}
\text{Min}  ~  f(x)& = (2\sqrt{2}x_1+x_2). l \\
\text{Subject to:} & \\
g_1(x)&=\dfrac{\sqrt{2}x_1+x_2}{\sqrt{2}x_1^{2}+2x_1x_2}P-\sigma \leq 0,\\
g_2(x)&= \dfrac{x_2}{\sqrt{2}x_1^{2}+2x_1x_2}P-\sigma \leq 0 \\
\end{split}
\end{equation}
$$
g_3(x)= \dfrac{1}{x_1+\sqrt{2}x_2}P-\sigma \leq 0, ~ \text{where}$$
$$ 0 \leq x_1, x_2  \leq 1, l=100cm, P=\sigma=2 KN/cm^{2}$$
\begin{figure}[ht]
	\centering
	\includegraphics[scale=0.75]{bt.eps}
	\caption{Structure of Three-bar truss}
\end{figure}

\begin{table}[!ht]
	\centering
	\caption{Comparison of optimal results for three-bar truss problem}
	\begin{tabular}{c|c|c|c|c}
		\hline
		Variables & PSO & QPSO & SP-QPSO & QPSO-CD  \\
		\hline
		$x_1$ &0.78911058  &0.788649  &  0.788796  & 0.788658\\
		$x_2$&  0.40702683 & 0.408322 &    0.407898& 0.40828488\\
		$g_1(x)$ & -6.6720e-06 & 1.6313e-07 & 6.4748e-06 &  9.00037e-06 \\
		$g_2(x)$ & -1.4655 & -1.4640 &   -1.4644 & -1.4640 \\
		$g_3(x)$ & -0.5345 &-0.5359  &    -0.5354 & -0.5359\\
		$f(x)$ & 263.89686 &  263.89584 &     263.89500 & 263.89465\\
		\hline
	\end{tabular}
\end{table}
The results are obtained by QPSO-CD are compared with its counterparts in Table 6. For three-bar truss problem, QPSO-CD is superior to optimal solutions  previously obtained in literature. The difference of best solution obtained by QPSO-CD among other algorithms is shown in Fig 9.

\begin{figure}[H]
	\centering
	\includegraphics[scale=0.55]{tbar.eps}
	\caption{Optimal results of PSO, QPSO, SP-QPSO and QPSO-CD algorithms for three bar truss problem}
\end{figure}

\subsection{Tension/Compression spring design problem}
The main aim is to lessen the volume \textit{V} of a spring
subjected to tension load constantly as shown in Fig 10. Using the symmetry of structure, there are practically three design variables ($x_1, x_2, x_3$), where $x_1$ is the wire diameter, the coil diameter is represented by $x_2$ and $x_3$ denotes the total number of active coils.    The mathematical formulation for this problem is described as:

\begin{equation}
\begin{split}
\text{Min}  ~  f(x)&= (x_3+2)x_2x_1^2, \\
\text{Subject to:} & \\
g_1(x)&=\dfrac{1-x_2^{3}x_3}{71785x_16{4}} \leq 0,\\
g_2(x)&= \dfrac{4x_2^{2}-x_1x_2}{12566(x_2x_1^{3}-x_1^{4})}+\dfrac{1}{5108x_1^{2}}-1 \leq 0,\\
g_3(x)&=  \dfrac{1-140.45x_1}{x_2^2x_3}\leq 0,\\
g_4(x)&=  \dfrac{x_2+x_1}{1.5}-1\leq 0,\\
\text{where} & \\
0.05 \leq x_1 & \leq 2,  0.25  \leq x_2, \leq 1.3, 2 \leq x_3  \leq 15, \\
\end{split}
\end{equation}

	\begin{figure}[H]
	\centering
	\includegraphics[scale=0.38]{ts.eps}
	\caption{Structure of tension/compression spring}
\end{figure}
\begin{table}[!ht]
	\centering
	\caption{Comparison of optimal results for tension spring design problem}
	\begin{tabular}{c|c|c|c|c}
		\hline
		Variables & PSO & QPSO & SP-QPSO & QPSO-CD  \\
		\hline
		$x_1$ &  0.0516 & 0.0524  & 0.05 & 0.0513   \\
		$x_2$&  0.3542 & 0.2505  &  0.25  & 0.2502  \\
		$x_3$ & 11.7942 & 2 &  2 & 2\\
		$g_1(x)$ & -2.3006e-02  & 0.93145095  & 0.93034756 & 4.11004e-06  \\
		$g_2(x)$ &  -5.6059e-03  &-0.17471558  & -0.16568318 & -0.17352479  \\
		$g_3(x)$ & -3.9057 & -50.67  & -48.180  &  -49.561 \\
		$g_4(x)$ & -0.7294 & -0.79986567 & -0.80   & -0.799 \\
		$f(x)$ & 0.01305 &0.00275  &  0.00250  &  0.00263\\
		\hline
	\end{tabular}
\end{table}

It has been observed that QPSO algorithm with Cauchy distribution and natural selection strategy is robust and obtains optimal solutions than
PSO and QPSO, shown in Table 7. The difference between best solutions found by QPSO-CD  ($f(x)=0.00263$) and other algorithms for tension spring design problem are reported in Fig 11.

\begin{figure}[H]
	\centering
	\includegraphics[scale=0.55]{tspg.eps}
	\caption{Results of PSO, QPSO, SP-QPSO and QPSO-CD methods for tension spring design problem}
\end{figure}

\subsection{Pressure vessel design problem}

Initially,  Kannan and Kramer \cite{kannan1994augmented} studied the pressure vessel design problem with the main aim to reduce the total fabricating cost. Pressure
vessels can be of any shape. For engineering purposes, a cylindrical design capped by hemispherical heads at both ends is widely used \cite{sandgren1990nonlinear}. Fig 12 describes the structure of pressure vessel design problem. It consists four design variables ($x_1, x_2, x_3, x_4$), where $x_1$ denotes the shell thickness $(T_s)$, $x_2$ is used for head thickness ($T_h$), $x_3$ denotes the inner radius (\textit{R}) and  $x_4$ represents the length of vessel (\textit{L}). The objective function and constraint equations are described as:
\begin{figure}[H]
	\centering
	\includegraphics[scale=0.62]{pv.eps}
	\caption{Design of pressure vessel}
\end{figure}

\begin{equation}
\begin{split}
\text{Min}  ~  f(x)& =  0.6224x_1x_3x_4 + 1.7781x_2x_3^{2}+  3.166x_1^{2}x_4 \\ & +19.84x_1^{2}x_3,  \\
\text{Subject to:} & \\
g_1(x)&= -x_1+0.0193x_3 \leq 0,\\
g_2(x)&=  -x_2+0.00954x_3\leq 0,\\
g_3(x)&= -\pi x_3^{2}x_4-\dfrac{4}{3}\pi x_3^{3}+1296000 \leq 0,\\
g_4(x)&= x_4-240 \leq 0,\\
\text{where} & \\
1\times 0.0625 & \leq x_1, x_2  \leq 99 \times 0.0625,  10  \leq x_3 ,x_4 \leq  200 \\
\end{split}
\end{equation}

\begin{table}[!ht]
	\centering
	\caption{Comparison of optimal results for Pressure vessel design problem}
	\begin{tabular}{c|c|c|c|c}
		\hline
		Variables & PSO & QPSO & SP-QPSO & QPSO-CD  \\
		\hline
		$x_1$ & 0.8125 &0.7783  & 0.7782   & 0.7776\\
		$x_2$&  0.4375& 0.3849 & 0.3845    & 0.3848\\
		$x_3$ & 42.0984 & 40.3289 & 40.3206 & 40.3278\\
		$x_4$ & 176.6365  & 199.8899 & 199.9988 & 199.8865 \\
		$g_1(x)$ & -4.500e-15 & 4.777e-05  & -1.242e-05  & 7.2654e-04\\
		$g_2(x)$ & -0.035880 & -1.62294e-04  & 1.58523e-04   & -7.2787e-05 \\
		$g_3(x)$ & -1.164e-10  & -97.39720071  & -63.63686942    & -0.734359\\
		$g_4(x)$ & -63.3634 & -40.1100 &  -40.0012  & -40.1135 \\
		$f(x)$ & 6059.714
		&5886.189  & 5885.268   & 5886.137\\
		\hline
	\end{tabular}
\end{table}

\begin{figure}[!ht]
	\centering
	\includegraphics[scale=0.55]{pvdp.eps}
	\caption{Optimal results of PSO, QPSO, SP-QPSO and QPSO-CD techniques for pressure vessel design problem}
\end{figure}

The optimal results of QPSO-CD is compared with the SP-QPSO, QPSO and PSO best results noted in the previous work, and are given in Table 8. The best solution obtained from QPSO-CD is  better than other algorithms as shown in Fig 13.

\section{Conclusion}
In this paper, a new hybrid quantum particle swarm optimization algorithm is proposed with natural selection method and Cauchy distribution. The performance of the proposed algorithm is experimented on four benchmark functions and the optimal results are compared with existing algorithms. Further, the QPSO-CD is applied to solve engineering design problems. The efficiency of QPSO-CD
is successfully presented with superiority  than preceding results
 for three engineering design problems: three-bar truss, tension/compression spring and pressure vessel. The efficiency of QPSO-CD algorithm is evaluated  by number of steps needed to reach the optimal region and proved that time complexity of proposed algorithm is lower in comparison to classical PSO.  In the context of convergence, the experimental outcomes shown that the QPSO-CD converge to get results closer to the superior solution.

\section*{Additional information}
\textbf{Competing interests:} The authors declare no competing interests.

\section*{Acknowledgement}
S.Z. acknowledges support in part from the National Natural Science
Foundation of China (Nos. 61602532),  the Natural Science Foundation of Guangdong Province of China (No. 2017A030313378), and the Science and Technology Program of Guangzhou City of China (No. 201707010194).