(r),\phantom{\Big(}\\
\end{array}
\right\}\quad r \gtrsim R_0,\end{aligned}$$ where $k_i(r) = \sqrt{2 \mu \left(E - U_{i}(r)\right)}/\hbar$ is the local wave vector, $U_{i}(r) = V_p(r) + \hbar^2 \ell(\ell+1)/(2 \mu r^2)$, and $\beta_i(r) = \int^r\!{\mathrm}{d}x\,k_i(x)$ is the WKB phase. Another possibility is to use the inhomogeneous Milne equation for parametrization [@Mies2000]. The total wave function at short range can be written as $$\label{Psi}
\mathbf{\Psi}(r,E)= \mathbf{A}(E) \left[\mathbf{\hat{f}}(r,E) + \mathbf{Y}(E,\ell) \mathbf{\hat{g}}(r,E) \right]$$ Here, $\mathbf{A}(E)$ is the amplitude, $\mathbf{\hat{f}}_{ij}=\hat{f}_i\delta_{ij}$, $\mathbf{\hat{g}}_{ij}=\hat{g}_i\delta_{ij}$ are diagonal matrices and $\mathbf{Y}(E,\ell)$ is the so-called quantum defect matrix, a crucial object in this method. At long range the solution of the problem can be expressed using energy normalized functions $$\begin{aligned}
\label{fg}
\left.
\begin{array}{ll}
f_i(r,E) \cong \sin \left(kr-\ell\pi/2+\xi_i\right)/\sqrt{k},\\
g_i(r,E) \cong \cos \left(kr-\ell\pi/2+\xi_i\right)/\sqrt{k},
\end{array}
\right\}\quad r \rightarrow \infty,\end{aligned}$$ where $\xi_i$ denotes the phase shift induced by the full potential $U_i$. The short- and long-range solutions can then be matched using the quantum defect functions $C(E,\ell)$ and $\tan\lambda(E,\ell)$ $$\begin{aligned}
\begin{array}{ll}
f_i(r,E) = & C_i^{-1}(E) \hat{f}_i(r)\\
g_i(r,E) = & C_i(E)(\hat{g}_i(r)+\tan \lambda_i(E) \hat{f}_i(r)).
\end{array}\end{aligned}$$ We notice that one can intuitively interpret $C$ and $\tan\lambda$ functions as a measure of deviation of the solution from WKB one. As a result it is clear that for energies high above threshold $C(E)\to 