We consider the non-symmetric (or Johnson-Nédélec) coupling to solve second order uniform elliptic partial differential equations defined in unbounded domains. We present a novel condition that ensures the ellipticity of the associated bilinear form, keeping track of its dependence on the linear combination coefficients of the interior domain equation with the boundary integral one. We show that an optimal ellipticity condition, relating the minimal eigenvalue of the diffusion matrix to the contraction constant of the shifted double-layer integral operator, is guaranteed by choosing a particular linear combination. For the numerical discretization of the scheme, we propose the coupling of a Virtual Element Method (VEM) with a Boundary Element Method (BEM), by using decoupled approximation orders. We provide optimal convergence error estimates, in the energy and in the weaker L2-norm, in which the VEM and BEM contributions to the error are separated. This allows taking advantage of the high order flexibility of the VEM to retrieve an accurate discrete solution by using a low order BEM.