r})=\frac{x+iy}{\sqrt{x^2+y^2+\xi_v^2}}
  =\frac{r e^{i\theta}}{\sqrt{r^2+\xi_v^2}}$$ with the core size $\xi_v=\xi$. To describe the system of two singly–quantized vortices positioned at ${\bf r}=\pm{\bf a}/2=\pm(a/2,0)$ we fit the gap function as follows: $$\begin{aligned}
\label{2VortMolecGap}
  \nonumber
  \Delta({\bf r})=\Delta_0f_1({\bf r}-{\bf a}/2)f_1({\bf r}+{\bf a}/2)\,.\end{aligned}$$

The image method allows us to reduce the problem of a vortex positioned at the point ${\bf r}=(a/2,0)$ near a flat boundary $x=0$ to the one describing a vortex-antivortex pair with the antivortex situated at ${\bf
r}=-(a/2,0)$. For the latter case we choose the gap function in the form: $$\begin{aligned}
\label{VaVVortMolecGap}
 \nonumber
 \Delta({\bf r})=\Delta_0f_1({\bf r}-{\bf a}/2)f_1^*({\bf r}+{\bf a}/2)\,.\end{aligned}$$

Certainly, the above simplest approximation for the gap function might not be enough to describe closely spaced vortex–vortex or vortex–antivortex pairs. However, the resulting LDOS patterns obtained below are mainly controlled only by the distance between the phase singularity points and appear to be weakly sensitive to the concrete profiles of the order parameter absolute value. The key features of the LDOS patterns appear to hold even in the limit of the zero $\xi_v$ value, i.e., even for phase vortices without any suppression of the gap value inside the cores.

Considering hereafter only the vortex configurations of the finite size comparable with several coherence lengths 