\section{Investigating patterns in the bias}
\label{appendix}
\begin{figure*}[htb!]
	\centering
	\includegraphics[width=0.5\textwidth]{figs/q05_inc90_mtot70_int.pdf}
	\includegraphics[width=0.45\textwidth]{figs/q05_inc90_mtot70_ext.pdf}
	\caption{The recovery of source parameters using Type-I{} (blue) and Type-II{} (red) templates for $q=2$, $M=70 M_\odot$, $\rho=50$, and $\iota=\pi/2$. The Type-II{} recoveries are consistent with the injection values, but the Type-I{} recoveries have correlated bimodalities $ t_c $, $ \delta $ and $ \iota $. One of these modes is consistent with the injected values, while the other mode is far away from the injected values. The distance posterior does not have a significant bias with respect to the injected value; however, the width of the posterior is larger than that of posterior in the Type-II{} case.}
	\label{fig:q05_inc90_mtot70}
\end{figure*}

\begin{figure*}[htb!]
	\centering
	\includegraphics[width=\textwidth]{figs/fdomain_amplitude_phase.pdf}
	\caption{Comparison of the waveform amplitude $ |h(f)| $ and phase $ \psi(f) $ (projected onto the LIGO-Hanford detector) for the maximum likelihood samples from the low-$ t_c $ and high-$ t_c $ modes. There is considerable agreement between the two waveforms.}
	\label{fig:fdomain_amplitude_phase}
\end{figure*}

For some cases in Fig. \ref{fig:SNR20_bias} and Fig. \ref{fig:SNR50_bias}, there are anomalous dips in the bias (quantified by $\mathcal{Q}_\mathrm{rec}$) at $\iota= \pi/2$. We investigate these below.

In order to illustrate the causes of these dips, we show the corner plot of source parameter recovery for a Type-II injection with $ q=2 $, $ M=70 M_\odot $, $ \rho=50 $ and $ \iota = \pi/2 $ with both Type-I and Type-II templates (Fig. \ref{fig:q05_inc90_mtot70}). The posteriors obtained with the Type-II template do not show any anomalous features, and the injection values lie near the peak of the posterior. However, this is not the case for the Type-I recovery. The first thing to note is that the $ t_c $ posterior is bimodal, thus also causing the $ \delta $ posterior to be bimodal since there is a correlation between the two parameters. This is perhaps not too surprising since, as we have shown earlier, Type-II effects distort the waveform away from GR considerably. One of the modes in the $ t_c - \delta $ plane includes the injected value, while the other is far away from it. This bimodality also effectively broadens the region where the posterior has significant support, thus causing the $ \mathcal{Q}_\mathrm{rec} $ to be shifted to lower values. This is the reason for the anomalous behaviour seen at $ \iota = \pi/2 $, and one can also see that a similar argument holds for the recovery of $ \iota $ as well. On the other hand, the $ d_L $ posterior peaks at the injected value, but the posterior width is larger than that  in the Type-I recovery case. This, again, causes a decrease in the value of $ \mathcal{Q}_\mathrm{rec} $ estimated for the $ d_L $ parameter. The points noted above show that $ \mathcal{Q}_\mathrm{rec} $ is not a good enough quantifier of the bias for these edge cases.

We have also checked that the likelihood values for samples at both the modes is very similar. To check whether the waveforms themselves at these two modes are similar, and to ensure sanity of the likelihood calculation and the sampling, we plot the frequency-domain amplitude and phase of the maximum likelihood waveforms (projected onto the LIGO-Hanford detector) from both the modes (Fig. \ref{fig:fdomain_amplitude_phase}). As one can clearly see, the waveforms match very well in the frequency domain. 

\section{Introduction}
	
	The sensitivities of the LIGO \& Virgo gravitational-wave (GW) detectors \cite{2015CQGra..32g4001L,2015CQGra..32b4001A} are constantly improving. So far, they have already confidently detected $90$ compact binary mergers events \cite{LIGOScientific:2021djp} and are expected to detect hundreds of such mergers in the upcoming observational runs \cite{KAGRA:2013rdx}. Accurate extraction of source parameters from these mergers is essential to interpret them in an astrophysical and cosmological context. Understanding and modeling physical effects due to the propagation of the wave between the source and the detector are extremely important in this regard. Ignoring these effects in the templates used to search and estimate parameters of a GW signal could cause problems in explaining or interpreting the underlying signals.
	
	One such propagation effect is gravitational lensing. GWs get lensed when the signals encounter mass inhomogeneities  during their journey to the detectors. If the intervening object is sufficiently massive, it can produce multiple images of the source through strong lensing \cite{Dai:2016igl,Ng:2017yiu}. The multiple images would differ in their magnifications and arrival times to the observer. The probability (quantified by the optical depth) for strong lensing is, however, small, e.g., $\sim 10^{-3} - 10^{-4}$ for galaxy lenses \cite{Wang:2021kzt}, mainly because it requires a strong alignment between the source, lens and the observer. Searches in GW data have not detected signatures of strong lensing so far \cite{Hannuksela:2019kle, LIGOScientific:2021izm}. Nevertheless, it has been anticipated that $\sim 1.3^{+0.6}_{-0.4} - 1.7^{+0.9}_{-0.6}$ binary black hole (BBH) lensed signals per year from galaxy lensing \cite{Wierda:2021upe} could be observed at the design sensitivities of advanced LIGO \& Virgo \cite{LIGOScientific:2014pky, VIRGO:2014yos}, where the uncertainty mainly arises from the BBH merger rate distribution. These prospects will only be enhanced with the addition of KAGRA \cite{Aso:2013eba, KAGRA:2020tym} and LIGO-India \cite{Saleem:2021iwi} to the detector network. Furthermore, third generation GW detectors like Cosmic Explorer (CE) \cite{Reitze:2019iox} and Einstein Telescope (ET) \cite{Sathyaprakash:2012jk} are expected to have an order of magnitude better sensitivity as compared to the current generation of detectors, thus potentially observing hundreds of thousands of mergers every year. This suggests an exciting time for doing science with strongly lensed GWs in future~\cite{Sereno:2011ty, Ding:2015uha, Baker:2016reh, Hannuksela:2020xor, Goyal:2020bkm, Xu:2021bfn}.
	
	In the geometric-optics limit, the lens equation picks solutions at extremal points of the time-delay (arrival time with respect to a reference time) surface over a family of trajectories. The images thus produced due to strong lensing can be of different types corresponding to the minima, saddle point and the maxima of the time-delay, called Type-I, Type-II and Type-III images, respectively. Type-I and Type-III images have positive parity (i.e., positive magnification) while Type-II images have negative parity. In addition, Type-II images receive a $\pi/2$ phase shift, while Type-I and Type-III images receive a phase shift of $0$ and $\pi$, respectively \cite{Dai:2017huk}.  Thus, Type-I and Type-III images can be mimicked by unlensed signals with an overall magnification. Type-II images, on the other hand, get distorted compared to their unlensed counterparts due to the non-trivial phase shift of $\pi/2$.  In the cases when the GW radiation is mainly composed of the dominant harmonics ($2, \pm 2$ modes) these distortions can be mimicked by Type-I or the unlensed GW signals by adjusting their parameters such as the coalescence phase and/or polarization angle \cite{Ezquiaga:2020gdt}. However, when the signal has significant contributions from higher harmonics, precession, or eccentricity, these distortions may not be reproduced by the Type-I signals \cite{Ezquiaga:2020gdt}, as also discussed in Sec.~\ref{sec:setup} here.  In such cases, the current strategies that LIGO-Virgo follows i.e., analyzing a detected signal under the hypothesis that the signal is unlensed, may lead to significant biases in the inferences of the parameters. The signals could also be missed by the standard LIGO-Virgo search pipelines.  
	
	Under the singular isothermal ellipse (SIE) approximation~\footnote{A more realistic approximation than the singular isothermal sphere (SIS).} of the galaxy lenses, it has been shown that the probability of seeing multiple images without Type-II image being one of them is less than $0.01\%$ \cite{Wang:2021kzt}. Hence, the optical depth for Type-II images is practically the same as that of strong lensing. Moreover, it was also shown that more than $90\%$ of sources with multiple images would have a Type-II image as the second brightest image \cite{Wang:2021kzt}. Thus, if a strongly lensed event pair is detected, there is a high chance that a Type-II image would be present.
	
	
	In this work, we investigate thoroughly the systematic effects that could arise when a Type-II (BBH) signal is recovered with the Type-I/unlensed templates. We first establish the cases where a single-event-based confident identification of the Type-II nature of a signal could be possible. This is done by choosing a threshold for the log Bayes factor between the two hypotheses $\textbf{A}$ and $\textbf{B}$, where  $\textbf{A}$ represents the hypothesis where we inject a Type-II signal and recover with Type-I templates, and $\textbf{B}$ represents the hypothesis when the Type-II injected signal is recovered with the Type-II template. We claim a Type-II signal is identified when the log Bayes factor $\ln \mathcal{B}$ $\geq 2$. We study the effect of strong lensing on a wide set of injections at total SNR of 20 and 50, using the noise-spectra of the LIGO-Virgo design sensitivities. We show that even at modest values of SNR, the Type-II nature of the image can be inferred for moderate to high values of the inclination. We also show that non-inclusion of the Morse phase in the parameter estimation templates for such signals could lead to large systematic effects, pushing the true values of the parameters beyond the $90\%$ credible intervals (CIs) of their posteriors. In other words, the parameters of the signals could be wrongly inferred if Type-II templates are not used during the parameter estimation (PE). This could, in turn, bias the astrophysical interpretation of the source, e.g., a lighter BH could be wrongly attributed to higher mass, thus causing difficulties in explaining their formation mechanisms. 
	
	One strategy to look for the strong lensing signatures in the data is by analyzing pairs of events with consistent sky locations and chirp masses \cite{Haris:2018vmn, Goyal:2021hxv}~\footnote{Ref. \cite{Caliskan:2022wbh} quantifies the false alarm probabilities associated with this strategy, and the difficulties that arise in identification of lensed events as a result.}. This is motivated by the belief that the strong lensing will not affect the frequency profile of the unlensed signal. This is indeed true for Type-I and Type-III images, and in some cases the Type-II images too if they are quasi-circular and contain only the dominant (22) mode. In general though, one would need to perform a joint analysis,  where  a pair of events are simultaneously analyzed  \cite{Liu:2020par, Lo:2021nae,Janquart:2021qov}. However, as we discussed before, in a realistic scenario if there is a lensed event pair in the data then it is highly likely that one of them would be a Type-II image. Then, given that we would not know \textit{a priori} the parameters of the underlying unlensed GW signal, we may not expect that there will be a consistency in the parameters (such as the chirp mass and the sky location parameters) extracted using the templates of the unlensed signals. Thus, in principle, one would need to do joint-PE search over the full data without any prior belief. This may become a difficult task for such techniques given their computational costs. On the other hand, doing a full search over data for just the individual Type-II images should be relatively easy and if detected, they can help narrow down the search for its corresponding pair since we would already know the parameters of the underlying unlensed signal. Moreover, in the cases where the lensed pair can not be detected due to detector down-time or data quality issues, signatures of lensing can be obtained from just one image and further used for downstream analyses such as estimating lensing rates.
	
	We note that a similar study has been performed in Ref. \cite{Wang:2021kzt}. However, the study therein is based on the analytical approximation of the Bayes factor, with emphasis on the identification of the Type-II signals. They also do not explicitly take into account the detector response of the GW signals, avoiding this by making an assumption that both polarizations are measurable. In our work, we perform a full parameter estimation study, emphasizing both the identification of the Type-II signals as well as the biases in the inference of the binary parameters. Thus, in this sense, our work builds upon their work along this particular goal. We also note that Ref. \cite{Janquart:2021nus} has performed a complementary analysis, showing that including higher-order modes in joint-PE allows identification of individual image types in a lensed pair.  
	
	This paper is organized as follows: in Sec.~\ref{sec:setup} we explain the effect of the strong lensing on GWs and describe the method followed in this work. In Sec.~\ref{sec:results}, we discuss the results based on a wide set of simulations tuned to our interests. Sec.~\ref{sec:conclusion} presents the conclusion and future work. 
	
	\begin{figure*}[htb!]
		\includegraphics[width=0.95\textwidth]{figs/Type2_lensed_waveform_Mtot100_q7.00_iota0.00.pdf}
		\includegraphics[width=0.95\textwidth]{figs/Type2_lensed_waveform_Mtot100_q7.00_iota1.05.pdf}
		\caption{The unlensed (GR), the Type-II, and the GR signal with a shift in the coalescence phase $\phi_c$ by $\pi/4$ for the parameters $M = 100 M_{\odot}$, $q=7$, $\iota=0$ (top panel) and $\iota=\pi/3$ (bottom panel). For the zero inclination (top panel), i.e., when mainly 22 mode contributes, the GR signal reproduces the Type-II image signal very well by adjusting its coalescence phase while for high inclination (bottom), Type-II signal is not reproduced by the GR signal fully, i.e., there are additional distortions left.  }
		\label{fig:lensed_waveforms}
	\end{figure*}
	
	\section{Set up}
	\label{sec:setup}
	\subsection{Theory}
	
	The problem of lensing of GWs has been looked at in-depth already \cite{PhysRevD.59.083001}.  In geometric optics limit (i.e., the short wavelength regime, valid for, e.g., lensing of GWs by galaxies), the effect of lensing on GWs can be derived from the Kirchhoff diffraction integral~\footnote{The Kirchhoff diffraction integral describes the generic lensing phenomenon.} by choosing points on the image plane corresponding to the extremal points of the time-delay \cite{1986ApJ...310..568B}. Since the extremal points can be either the minimum, saddle or the maximum, we can have three types of images, namely, Type-I, Type-II, Type-III, respectively. In this limit, the Kirchhoff integral reduces to the Gaussian integral around these extremal points and the result is the following amplification factor for $j^{\rm{th}}$ image,
	
	\begin{equation}
	F_{j}(f) = \sqrt{\mu_j}  \exp[\mathrm{i} 2\pi f \delta t_j-\text{sign}(f) n_j\pi]
	\end{equation}
	where $f$ is the observed GW frequency, $\mu_j$ and $\delta t_j$ represent the magnification and time delay of the $j^{\rm{th}}$ image, and the extra phase shift $\Delta \phi_{j}= \text{sign}(f) n_j\pi$ results from the integration of the complex Gaussian function, where $n_j$, called the Morse index, takes the values $0$, $1/2$ or $1$ for Type-I, II or III images, respectively. The lensed GW signal is then given by convolving the unlensed signal with the amplification factor,  
	
	\begin{equation}
	\Tilde{h}_{+,\times,\, j}^{L} (f) = F_{j}(f) \ \Tilde{h}_{+,\times}(f), 
	\label{eq:amp_fact}
	\end{equation}
	As we can see, the Type-I image just receives a linear phase compared to the unlensed signal. This linear phase is practically un-important for all of our purposes since this can not be measured by the PE. The Type-III image  receives an additional overall negative sign but otherwise is again the same as the unlensed signal. The Type-II image, on the other hand,  depending on the sign of the frequencies can receive $+\pi/2$ or $-\pi/2$ phase shifts. As a consequence, the time-domain lensed signal is given by the Hilbert transformation of the Eq.~\eqref{eq:amp_fact} rather than the simple inverse Fourier transformation of the unlensed signal,
	
	\begin{equation}
	h_{+,\times}^{L,\mathrm{II}} (t) = - \int_{-\infty}^{\infty} df \mathrm{sign} (f) i e^{-i2\pi f t}  \Tilde{h}_{+,\times}(f)
	\label{eq:time_dom_lensed_wav}
	\end{equation}
	
	For our purpose in this work, we directly use the Eq.~\eqref{eq:amp_fact} for the lensed template construction because the GW data analysis is done in the frequency domain itself. Nevertheless, Eq.~\eqref{eq:time_dom_lensed_wav} serves a better purpose for the visual inspections of the lensing effects and thus we use it just for the sake of demonstrations. Also, the standard GW data analysis employs only the positive frequencies because, for a real signal, the negative frequencies do not provide any extra information and hence only $+ \pi/2$ phase shift is relevant for us. We now discuss below the implications of the $\pi/2$ phase shift for the Type-II images when they have different structures. 
	
	A general GW can have different harmonic contents. One way to visualize them is via spherical mode decomposition,
	\begin{equation}
	h_{+}(t) -i h_{\times}(t) = \sum_{\ell \geq 2}\sum_{m=-\ell}^{\ell} h_{\ell m}(t; \bm{\lambda}) _{-2}Y_{\ell m}(\iota, \phi_{c})
	\label{eq:mode_decomp}
	\end{equation}
	where $_{-2}Y_{\ell m}$ \cite{Pan:2011gk} denotes the spin-2 weighted spherical harmonics which are function of the inclination angle $\iota$~\footnote{The angle between the total angular momentum $\vec{\boldsymbol{J}}$ of the binary and the observer.} and the coalescence phase $\phi_c$. The spherical harmonics basically separate the angular part of the GW radiation from its radial part and thus the modes $h_{\ell m}(t)$ depend just on the intrinsic parameters ($\bm{\lambda}$) of the binary, e.g., component masses ($m_{1,2}$) and spins ($\chi_{1,2}$) if the binary is circular. For non-precessing BBH binaries in circular motion, the $m<0$ modes are related to the $m>0$ due to the reflection symmetry of the GW radiation about the binary plane,
	\begin{equation}
	h_{\ell, -m} = (-1)^{\ell} h^{*}_{\ell m}
	\end{equation}
	This helps simplify Eq.~\eqref{eq:mode_decomp} in terms of only the $m>0$ modes and thus throughout this paper we use $m>0$ modes to denote the full mode contents. On a bit further simplification of Eq.~\eqref{eq:mode_decomp} writing $h_{\ell m}= A_{\ell m}e^{i \phi_{\ell m}}$, one would obtain (see Appendix C of \cite{Mehta:2017jpq}),
	
	\begin{align}
	h_{+}(t) &=  \sum_{\ell \geq 2}\sum_{m\geq 0} f_{\ell m}^{+}(\iota) A_{\ell m}(t) \cos[\phi_{\ell m}(t) + m \phi_c] \label{eq:plus}\\
	h_{\times}(t) &=  \sum_{\ell \geq 2}\sum_{m\geq 0} f_{\ell m}^{\times}(\iota) A_{\ell m}(t) \sin[\phi_{\ell m}(t) + m \phi_c]
	\label{eq:cross}
	\end{align}
	where $f_{\ell m}^{+}(\iota)$ and $f_{\ell m}^{\times}(\iota)$ are functions of the inclination angle encoding the magnitude of the spherical harmonics $_{-2}Y_{\ell m} (\iota)$. Their explicit expressions are not relevant for the discussion here. The relative contribution of the modes $\ell m$ to the polarizations $h_{+,\times}$ depend on the specific choice of the parameter $\bm{\lambda}$ and the inclination angle ($\iota$). For example, for non-precessing BBH binaries with even the moderate mass ratio $q=m_1/m_2 \geq 4$ the contribution of higher order modes (HMs) becomes significant, neglecting which can have consequences for detection and parameter estimation, e.g., the systematic bias, loss in the detection volume, etc.~\cite{PhysRevD.93.084019,PhysRevD.90.124004, Capano:2013raa}. 
	
	When a GW signal reaches the detectors, the strain induced on a particular detector is given by,
	
	\begin{equation}
	h = F_{+}(\theta, \phi, \psi) h_{+}(t) + F_{\times}(\theta, \phi, \psi) h_{\times}(t)
	\label{eq:strain}
	\end{equation}
	where $F_{+}$ and $F_{\times}$ are called the antenna pattern functions which arise as a result of the detector's response to the GW signal. They represent the angular sensitivity of that detector and hence is a function of the sky angles $\theta$ and $\phi$; these parameters denote the location of the binary in the sky. The angle $\psi$ denotes the freedom in the choice of fixing the binary plane with respect to the detector plane and is called the polarization angle. The explicit expressions of the antenna pattern functions for the LIGO \& Virgo detectors are,
	
	\begin{align}
	F_{+} = &\dfrac{1}{2} \big[1 + \cos^2(\theta)\big ]\cos(2\phi)\cos(2\psi) \nonumber \\
	&-\cos(\theta)\sin(2\phi)\sin(2\psi),  \\
	F_{\times} = & \dfrac{1}{2}\big[1 + \cos^2(\theta)\big ] \cos(2\phi)\cos(2\psi) \nonumber \\
	&+\cos(\theta)\sin(2\phi)\cos(2\psi) 
	\end{align}
	
	On combining the Eqs.~\eqref{eq:plus}-\eqref{eq:cross} and the Eq. \eqref{eq:strain} one would obtain, 
	
	\begin{equation}
	h = \sum_{\ell \geq 2}\sum_{m\geq 0} \mathcal{A}_{\ell m} (t) \cos[\phi_{\ell m}(t) + m\phi_c - \chi_{\ell m}]
	\label{eq:final_strain}
	\end{equation}
	
	where,
	
	\begin{align}
	\chi_{\ell m} = \tan^{-1}\Bigg( \dfrac{F_{\times}(\theta, \phi, \psi)f^{\times}_{\ell m}(\iota)}{F_{+}(\theta, \phi, \psi) f^{+}_{\ell m}(\iota)}\Bigg)
	\end{align}
	
	\begin{equation}
	\mathcal{A}_{\ell m} (t) = A_{\ell m}(t) |F_{+}|[1 + \tan^2(\chi_{\ell m})]^{1/2} 
	\end{equation}
	Eq.~\eqref{eq:final_strain} is just a bit more simplified version of the Eq.~\eqref{eq:strain}. Now by taking the Fourier transform of the above equation and then applying the Hilbert transformation in the Eq.~\eqref{eq:time_dom_lensed_wav} one can easily show that the strains produced by the lensed signals in the time-domain are given as follows: 
	
	For Type-I images, there is only an additional magnification factor $\mu_{I}$, i.e.,
	
	\begin{equation}
	h_{\rm{I}} = \sum_{\ell \geq 2}\sum_{m\geq 0} |\mu_{\rm{I}}|^{1/2} \mathcal{A}_{\ell m} (t_{\rm{I}}) \cos[\phi_{\ell m}(t_{\rm{I}}) + m\phi_c - \chi_{\ell m}]
	\label{eq:final_strain_type1}
	\end{equation}
	
	For Type-II images there is an additional $\pi/2$ phase shift,
	
	\begin{equation}
	h_{\rm{II}} = \sum_{\ell \geq 2}\sum_{m\geq 0} |\mu_{\rm{II}}|^{1/2} \mathcal{A}_{\ell m} (t_{\rm{II}}) \cos\Big[\phi_{\ell m}(t_{\rm{II}}) + m\phi_c - \chi_{\ell m} + \pi/2\Big]
	\label{eq:final_strain_type2}
	\end{equation}
	
	while for Type-III images,
	
	\begin{equation}
	h_{\rm{III}} = \sum_{\ell \geq 2}\sum_{m\geq 0} |\mu_{\rm{III}}|^{1/2} \mathcal{A}_{\ell m} (t_{\rm{III}}) \cos[\phi_{\ell m}(t_{\rm{III}}) + m\phi_c - \chi_{\ell m} + \pi]
	\label{eq:final_strain_type3}
	\end{equation}
	where $t_{i} = t + \delta t_{i}$ for the $i^{\rm{th}}$ image type.
	As expected, the Type-I and Type-III induced strains are just the rescaled versions of the ones caused by the unlensed GW signals. The extra $\pi$ phase shift in the Type-III strains which causes just a global sign flip can be easily mimicked by the polarization angle $\psi$ through $\chi_{\ell m}$ term under the change $\psi \rightarrow \psi + \pi/2$. Thus, for Type-III image signals, the PE with unlensed GW templates would yield a bias in the polarization angle by $\pi/2$. 
	
	\begin{figure*}[t]
		\centering
		\includegraphics[width=\textwidth]{figs/lnB.pdf}
		\caption{BBH injections with a total SNRs $\rho= 20$ (upper row) and $\rho=50$ (lower row) shown in the $\iota$-$q$ plane for different total masses. The color bar shows the log Bayes factor ($\ln \mathcal{B}$) between the Type-II and the Type-I recoveries. Injections which cross our threshold of $\ln \mathcal{B} \geq 2$ are plotted using triangles (i.e., distinguishable image types) while the others are plotted using circles.}\label{fig:Bayes_Factor}
	\end{figure*}
	
	The Type-II induced strain is interesting since  the $\pi/2$ phase shift can not be easily absorbed into other known parameters unless the signal is quasi-circular, non-precessing and comprised of just the dominant $\ell m =22$ mode, in which case, the coalescence angle $\phi_c$ can easily absorb it by adjusting itself, $\phi_c \rightarrow \phi_c + \pi/4$, as should be clear from Eq.~\eqref{eq:final_strain_type2} when $m=2$. One can note, however, that the polarization angle $\psi$ can also try to adjust itself such that $\Delta \chi_{\ell m} = -\pi/2$ and thus it is not always the case that only $\phi_c$ will get biased during the PE as we will show later. When HMs also start contributing to the strain, the different $\ell m$ modes would require different shifts in the $\phi_c$, i.e.,
	
	\begin{equation}
	\phi_c \rightarrow \phi_c + \dfrac{\pi}{2m}
	\end{equation}
	in order to mimic the $\pi/2$ lensing phase shift. Thus, the strain produced by Type-II image signals would look distorted compared to its unlensed counterpart. Fig.~\ref{fig:lensed_waveforms} (bottom panel) shows an illustrative example of this case with the binary parameters $M=100 M_{\odot}$, $q=7$ and $\iota = \pi/3$. Since the inclination and the mass ratio are high here, HMs also contribute significantly and thus there are additional distortions that are not mimicked by the change in the coalescence phase. In the top panel, however, we also show the case with the zero inclination for the same binary configuration. For very small inclinations, the GW signals mainly contain the 22 mode and thus we can see that the Type-II image is very well mimicked by the GR signal with the transformation $\phi_c \rightarrow \phi_c + \pi/4$, as expected.
	
	We now explore how useful these distortions in the Type-II image strains can be in separating such signals out from a catalog of detected BBH events, and at the same time  how they would affect the inference of parameters if the correct template is not used. Below we discuss the details of the simulations that we use to address these questions.
	
	\begin{figure*}[htb!]
		\centering
		\includegraphics[width=0.5\textwidth]{figs/q05_inc30_mtot100_int.pdf}
		\includegraphics[width=0.45\textwidth]{figs/q05_inc30_mtot100_ext.pdf}
		
		\caption{The recovery of source parameters using Type-I{} (blue) and Type-II{} (red) templates for $q=2$, $M=100 M_\odot$, $\rho=20$, and $\iota=\pi/6$. The plots show recovery of $\mathcal{M}$,  $q$, luminosity distance $d_L$, inclination $\iota$, sky locations $\delta$ (dec) and $\alpha$ (RA), and time of arrival $t_c$ in the geocenter frame. The orange lines show the injected value. There is no appreciable difference between the posterior distributions for both the recoveries due to the low inclination and less SNR in the higher modes.}
		\label{fig:q05_inc30_mtot100}
	\end{figure*}
	
	\begin{figure*}[htb!]
		\centering
		\includegraphics[width=0.5\textwidth]{figs/q05_inc75_mtot100_int.pdf}
		\includegraphics[width=0.45\textwidth]{figs/q05_inc75_mtot100_ext.pdf}
		
		\caption{The recovery of source parameters using Type-I{} (blue) and Type-II{} (red) templates for $q=2$, $M=100 M_\odot$, $\rho=20$, and $\iota=5 \pi/12$. The axis labels are the same as described in Fig. \ref{fig:q05_inc30_mtot100}. There is now a bias in the inclination and the distance, and the injected value is not recovered within 1$\sigma$. The bias in the distance also biases the recovery of the source frame masses. There are also some insignificant biases in the recovery of the sky locations and time of arrival.}
		\label{fig:q05_inc75_mtot100}
	\end{figure*}
	
	
	\begin{figure*}[t]
		\centering
		\includegraphics[width=\textwidth]{figs/bias_snr20_rec.pdf}
		\caption{ Each panel shows the quantiles ($\mathcal{Q}_{\mathrm{rec}}$) of the injection value recoveries against the inclination angle ($\iota$) of the injected BBH systems for a given total mass ($M$) (labeled on the right side of the panel) and SNR $\rho=20$. The injections are the Type-II lensed non-spinning BBH signals and the recoveries are with Type-I template (solid lines) and the Type-II template (dashed lines) accounting for the additional phase shift of $\pi/4$. At the top of the first panel, we label the parameters of the BBH systems that we are interested in. In each plot the different colors denote the different mass ratios of the BBH systems. A lower quantile value denotes a more accurate inference of a parameter. }\label{fig:SNR20_bias}
	\end{figure*}
	
	\begin{figure*}[t]
		\centering
		\includegraphics[width=\textwidth]{figs/bias_snr50_rec.pdf}
		\caption{Same as Fig.~\ref{fig:SNR20_bias} but at $\rho = 50$.} \label{fig:SNR50_bias}
	\end{figure*}
	
	\subsection{Simulations and Bayesian inference}
	Measurement of parameters from a GW signal amounts to estimating the posterior probability density function (PDF) $ p (\theta| d) $, where $ \theta $ is the set of parameters that describes a GW signal, and $ d $ is the data stream from the interferometers in the network. According to Bayes theorem, the posterior PDF is given by,
	\begin{equation}\label{key}
	p(\theta | d) = \dfrac{\mathcal{L}(d | \theta) \ \pi(\theta)}{\mathcal{Z}} \qq{.}
	\end{equation}
	Here, $ \mathcal{L}(d | \theta) $ is the likelihood function of the data given the parameters, $ \pi(\theta) $ is the prior distribution for the parameters, and the normalization factor $ Z $ is the evidence, given by 
	\begin{equation}\label{key}
	Z = \int_{\theta} \dd{\theta'}\mathcal{L}(d | \theta') \ \pi(\theta') \qq{.}
	\end{equation}
	The evidence can be thought of as the likelihood function marginalized over all the parameters. The evidence is not a useful quantity by itself but is very useful while comparing different models/hypotheses of describing the data. In the context of this work, the two hypotheses that we consider are 
	\begin{align*}\label{key}
	&\text{A: the signal is unlensed} \\
	&\text{B: the signal is Type-II lensed}. 
	\end{align*}
	In order to understand which hypothesis fits the data better, we formulate the Bayes factor $ \mathcal{B}^{\text{Type-II}}_{\text{unlensed}} $, which is defined as the ratio of the evidences of the two hypotheses
	\begin{align}
	\mathcal{B}^{\text{Type-II}}_{\text{unlensed}} &= \dfrac{Z_\text{Type-II}}{Z_\text{unlensed}}\\
	\implies \ln \mathcal{B} &= \ln \mathcal{B}^{\text{Type-II}}_{\text{unlensed}} = \ln Z_\text{Type-II} - \ln Z_\text{unlensed}
	\end{align}
	We simulate Type-II image strains using Eq.~\eqref{eq:amp_fact} for total mass $M=m_1+m_2$ $\in \{30, 70, 100, 140\} \,{\rm M_\odot}$, mass ratio $q = m_1/m_2 \in \{2, 3, 5, 7, 10\}$, and inclination $i \in \{\pi/6, ~\pi/4, ~\pi/3, ~5\pi/12, ~\pi/2\}$ of the non-spinning BBHs. 
	We scale the luminosity distance of these injections to total signal-to-noise ratios (SNRs), $ \rho $, of 20 and 50 in a network with the LIGO-Hanford and LIGO-Livingston detectors at the Advanced LIGO sensitivity, and the Virgo detector at the Advanced Virgo sensitivity. We choose a realization of the detector noise that is exactly zero at all frequencies (called the ``zero-noise'' realization); this can be thought of as the most probable realization of the noise. This choice allows us to isolate the effects of waveform systematics on parameter estimation from the noise systematics. However, a fully-rigorous analysis would require adding a nonzero noise realization to these injections in order to make it directly relevant for the LIGO-Virgo observations. We leave this for future studies. Nevertheless, the SNRs (such as 50) considered in our work are expected to be high enough to minimize the effect of noise systematics in our results.
	
	We perform full bayesian parameter estimation on these injections (simulations) using Type-I (Eq.~\eqref{eq:final_strain_type1}) and Type-II (Eq.~\eqref{eq:final_strain_type2}) templates, keeping the priors on the source parameters equivalent between the recoveries for each injection. The priors we choose are uniform in the detector-frame chirp mass~\footnote{The chirp mass $\mathcal{M}$ is a mass parameter which describes the inspiral part of the BBH signals at leading order	.}  and mass ratio of the binaries, uniform in comoving volume for the luminosity distance, and uniform in sky location and inclination of the binaries. The parameter estimation runs are performed using \textsc{bilby} and \textsc{bilby\_pipe} \cite{Ashton:2018jfp,Romero-Shaw:2020owr} coupled with the dynamic nested sampler \textsc{dynesty} \cite{Speagle:2019ivv}. We use \textsc{IMRPhenomXHM} \cite{Pratten:2020ceb} as our waveform approximant for the runs; this approximant contains the 33, 44, 21 and 32 modes of the gravitational wave signal in addition to the dominant 22 mode. 
	
	\section{Results}
	\label{sec:results}
	
	Fig.~\ref{fig:Bayes_Factor} shows the Bayes factor between the Type-II and Type-I recoveries of Type-II injected BBH signals at the $\rho=20$ and 50, respectively, for various total masses. The x-axis denotes the mass ratio while the y-axis denotes the inclination angle. We choose a threshold of $\ln \mathcal{B}=2$ for the  distinguishability of the Type-II recovery against the Type-I recovery of Type-II injected signals \cite{Jeffreys:1939xee, Kass:1995loi}. To zeroth order, a higher threshold in the Bayes factor would require the lensed events to be observed at higher SNRs keeping other parameters fixed. The triangle shapes in  Fig.~\ref{fig:Bayes_Factor} denote the signals which pass this threshold. We can see that at the $\rho=20$, independently of the total mass and the mass ratio, when the inclination $\iota \gtrsim 5\pi/12$ the Type-II image signals become distinguishable from the Type-I (or the unlensed signals). The same holds at the $\rho=50$ but now the threshold in the inclination decreases to $\sim \pi/3$ as one would have expected since at higher SNRs the effect of HMs is larger. Since $\ln \mathcal{B} \sim \rho^2$, we expect that at even higher SNRs (e.g., 100 or larger), even smaller inclinations ($\sim \pi/4$) should allow us to distinguish the Type-II images. This means that third-generation detector networks, where SNRs of $\sim 100$ could be common, would easily allow us to distinguish the Type-II images even when the inclination is low. This prediction is also consistent with results from Ref. \cite{Wang:2021kzt}.  
	
	
	For the distinguishable Type-II image signals, we then explore the effect on the recovery of  source parameters if they are recovered with Type-I image (or equivalently, unlensed) templates. Fig.~\ref{fig:q05_inc30_mtot100} compares the recovery of an injected Type-II image having $\rho=20$, $ M=100 M_{\odot} $, $ q=2 $, and $ \iota=\pi/6 $ with both Type-I and Type-II templates. We see that, for both recoveries, the injected values lie well within the posterior PDFs and there is no discernible differences between the two recoveries. However, increasing $ \iota $ to $ 5 \pi / 12 $ (Fig.~\ref{fig:q05_inc75_mtot100}) keeping the other parameters fixed causes the injected masses, distance and inclination values to lie in a region with low posterior probability if the recovery is done with a Type-I template. This shows that the distortions in the waveform can cause significant biases in parameter estimation in certain regions of the parameter space. We investigate this further in Fig.~\ref{fig:SNR20_bias} and \ref{fig:SNR50_bias}, where we show results of recoveries with Type-I image templates (Eq.~\eqref{eq:final_strain_type1}) denoted by solid lines and with Type-II image templates (Eq.~\eqref{eq:final_strain_type2}) denoted by the dashed lines at SNR of 20 and 50, respectively, for different total masses of the Type-II injected signals. In each plot, the x-axis denotes the inclination and different colors represent the different mass ratios. On the y-axis we show the quantiles \ensuremath{\mathcal{Q}_\mathrm{rec}}{} of the recovery of the injected values for the significantly affected parameters. For example, a $\ensuremath{\mathcal{Q}_\mathrm{rec}}=90\%$ denotes that the injected value of the parameter is recovered at a value that forms one edge of the $ 90\% $ posterior area; in other words, the injection value is recovered in the tail of the posterior. Such a case would represent a Type-II signal which would provide a biased estimate of this parameter upon recovery with the unlensed template. 
	
	We can see from the Fig.~\ref{fig:SNR20_bias} that, at $\rho=20$, for the Type-II signals with lower total mass ($M=30M_{\odot}$), the intrinsic parameters chirp mass $\mathcal{M}$ and mass ratio $q$ are recovered at nearly the same quantiles, i.e., the phase shift does not affect their recoveries. This could be because, as explained before, the HMs are not very relevant for the inspiral-domainted signals at low SNRs.  At higher SNRs, however, the difference between the quantiles of their recoveries becomes significant at high inclination angles ($\iota\gtrsim 5 \pi/12$) (Fig.~\ref{fig:SNR50_bias}). At such high SNRs, for higher total mass systems, the differences become significant even at the inclination $\iota\sim \pi/3$. This threshold in the inclination would decrease with increasing SNR. 
	We thus expect that the chirp mass and mass ratio would be affected significantly for the 3G detectors when $\rho \sim 100$ or greater could be achievable. We note that even at $\rho=50$, there are multiple binaries for which the quantiles of the chirp mass and mass ratio recoveries exceed 90\% and thus their estimates are completely biased (see Fig.~\ref{fig:SNR50_bias}).
	
	The inferences of the distance ($d_L$) and the inclination also get affected significantly at higher SNRs. At lower SNR (Fig.~\ref{fig:SNR20_bias}), there is not much difference between their recovery quantiles, however, the quantiles do get pushed beyond 90\% for a range inclination angles ($\pi/3 \lesssim \iota \lesssim 5 \pi/12$). At higher SNR (Fig.~\ref{fig:SNR50_bias}), on the other hand, the binaries even with the inclination angle of $\sim \pi/4$ face significant differences. This suggests that a given Type-II lensed event in the 3G detectors era will have biased estimates of its distance and inclination too unless the unlensed signal is really from a nearly face-on binary.
	
	
	We also find the declination ($\delta$) and the coalescence time ($t_c$) being affected very much, again at higher SNR. We can see from Fig.~\ref{fig:SNR50_bias} that higher total mass systems with the  angle of $\sim \pi/3$ onward have completely biased estimates of  these parameters. This threshold in the inclination would again decrease at even higher SNRs and we again expect that almost all of the Type-II lensed events in the 3G era and beyond will have this problem if recovered with the unlensed templates. In other words, the sky location inferred from these events will be highly inaccurate. 
	
	There are some seemingly anomalous features at $ \iota=\pi/2 $ in Fig. \ref{fig:SNR20_bias} and Fig. \ref{fig:SNR50_bias}, where the value of $  \mathcal{Q}_\mathrm{rec}$ drops as compared to $ \iota = 5 \pi/12 $. These occur due to bimodalities in $ t_c $, $ \delta $ and $ \iota $ posteriors (see Appendix \ref{appendix} for more details), thus rendering $ \mathcal{Q}_\mathrm{rec}$ to be an insufficient quantifier for the bias at $ \iota = \pi/2 $. 
	
	\section{Conclusion and Future Work}
	\label{sec:conclusion}
	
	Strong lensing produces three types of images; Type-I, Type-II and Type-III corresponding to the minima, saddle and maxima of the total arrival time of the lensed GWs. While Type-I and III images do not affect the frequency profiles of the unlensed BBH signals the Type-II image can, depending on the morphology of the signals, e.g., if the unlensed BBH signals have contributions from HMs or have support for the precession or eccentricity or any combination of these; for such signals the Type-II image strains receive additional distortions due to the phase shift ($\pi/2$) caused by the strong lensing. In this work we then explored the applicability of these distortions.
	
	Using a set of full parameter estimation recoveries for Type-II{} lensed injections on a varied parameter space, we showed that it would be possible to ascertain the Type-II{} nature of the most lensed BBH signals at high SNRs ($\gtrsim 100$). This becomes very relevant for the 3G detectors and beyond where such high SNR events could be frequently expected. This study implies that we should be able to tell if an event is Type-II{} lensed or not just by the observation of a single image. 
	
	We then showed that the identifiable Type-II images, which are of our interest here, will have sufficient distortions caused by the lensing phase shift $\pi/2$ that when they are recovered with Type-I image templates we would recover significantly biased parameters of the underlying BBH signals. We quantified these biases for a range of inclinations, total mass, and mass ratios of the binary.  Our results suggest that the recovery of parameters such as the chirp mass, mass rattio, the coalescence time, luminosity distance, inclination and the sky position of the binary becomes significantly biased and in many cases completely biased. Since the first-cut search for lensed pairs involves demanding consistency between the sky positions of the images, we suggest that parameter estimation should be strictly performed using Type-II templates on all signals, and the overlaps should be calculated between both Type-I{} and Type-II{} recoveries for both signals in the pair.
	
	In this work, we have considered the nonspinning quasicircular BBH signals in order to concentrate on the effects of only including higher harmonics of the radiation. It would be natural to extend this study to generic spinning BBH signals with eccentricity. We, however, expect that precession and/or eccentricity would cause similar biases. We have already verified this from a few precessing BBH Type-II injections but a thorough study will shed more light on this. We leave these investigations for future work.
	
	
	\section{Acknowledgements}
	We thank Jose Maria Ezquiaga, Justin Janquart and Otto Hannuksela for discussions. We also thank Bala Iyer, Parameswaran Ajith and other members of the Astrophysical Relativity group at ICTS for feedback, and Sanskriti Chitransh for a careful reading of the draft. Computations were performed on the Alice cluster at ICTS. AV's research is supported by the Department of Atomic Energy, Government of India, under Project No. RTI4001. This work makes use of \textsc{NumPy} \cite{vanderWalt:2011bqk}, \textsc{SciPy} \cite{Virtanen:2019joe}, \textsc{Matplotlib} \cite{Hunter:2007}, \textsc{jupyter} \cite{jupyter}, \textsc{dynesty} \cite{Speagle:2019ivv}, \textsc{bilby} \cite{Ashton:2018jfp}, \textsc{PyCBC} \cite{pycbc}, \textsc{LALSuite} \cite{lalsuite} and \textsc{PESummary} \cite{Hoy:2020vys} software packages.

    