Inspired by the works of Zagier (Eisenstein series and the Riemann zeta function, Automorphic forms, representation theory and arithmetic, Tata Institute of Fundamental Research, Bombay 1979, Springer, New York, 1981) and Sarnak (Commun Pure Appl Math, 34:719–739, 1981), we study the probability measures $$\nu (t)$$ with support on the flat tori which are the compact orbits of the maximal unipotent subgroup acting holomorphically on the positive orthonormal frame bundle $$\mathcal {F}({M}_D)$$ of 3-dimensional hyperbolic Bianchi orbifolds $${M}_D=\mathbb {H}^3/\widetilde{\Gamma }_D$$, of finite volume and with only one cusp. Here $$\widetilde{\Gamma }_D\subset \mathbf {PSL}(2,\mathbb {C})$$ is the Bianchi group corresponding to the imaginary quadratic field $$\mathbb {Q}(\sqrt{-D})$$. Thus $$\widetilde{\Gamma }_D$$ consist of Mobius transformations with coefficients in the ring of integers of $$\mathbb {Q}(\sqrt{-D})$$. If $$l \in \mathbb {N}$$, $$k,m \in \mathbb {Z}$$ are such that $$k,m \in [-l,l]$$, the appropriate Eisenstein series $$\widehat{E}_{km}^l ( g,s)$$ (which are defined and analytic for $$\mathrm {Re}(s) \ge 1$$) admit an analytic continuation to all of $$\mathbb {C}$$, except when $$l=k=m=0$$ in which case there is a pole for $$s=1$$. Using this fact we show that the elliptic curves which are the compact orbits of the complex horocycle flow $$h_T:\mathcal {F}({M}_D)\longrightarrow \mathcal {F}({M}_D)$$ ($$T \in \mathbb {C}$$) are expanded by the real geodesic flow $$g_t, t\in \mathbb {R}$$, and they become equidistributed in $$\mathcal {F}({M}_D)$$ with respect to the normalized Haar measure as $$t \longrightarrow \infty $$. This follows from the equidistribution of Lebesgue probability measures on compact leaves of the horocycle foliations in the orthonormal frame bundle of $$M_D$$, which is equal to the quotient $$\mathbf {PSL}(2,\mathbb {C})/\widetilde{\Gamma }_D$$. The same equidistribution property occurs for the spin bundle of $$M_D$$ which is the homogeneous space $$\mathbf {SL}(2,\mathbb {C})/\Gamma _D$$, where $$\Gamma _D$$ is the Bianchi subgroup in $$\mathbf {SL}(2,\mathbb {C})$$ consisting of matrices with entries in the ring of integers of $$\mathbb {Q}(\sqrt{-D})$$. Our method uses the theory of spherical harmonics in the unit tangent bundle orbifold $$T_1(M_D)={\mathbf {SO}}(2){\backslash }\mathbf {PSL}(2,\mathbb {C})/\widetilde{\Gamma }_D=T_1(\mathbb {H}^3)/{\Gamma ^*_D}$$, where $${\Gamma ^*_D}$$ is the action of $${\widetilde{\Gamma }}_D$$ on the unit tangent bundle of $$\mathbb {H}^3$$ via the differential of the elements of $$\widetilde{\Gamma }_D$$.