It is embarrassing to write my own answer to my question, but I do so for future readers seeking for a solution of similar integrals.
The integral $I(a,\theta)$, defined in the question turns out to be written in terms of error functions and Owen's T function, (i.e., it implies that the integral has something to do with bivariate normal distribution).
The answer is
$$
I(a,\theta) = 2 \pi T\left(\sqrt{2} a \cot (\theta ),\tan (\theta )\right)+\frac{1}{2}
\pi \text{erf}(a) (\text{erf}(a \cot (\theta ))-1),
$$
where $T$ is Owen's T function and $a$ is assumed to be positive.
Edit:
Let me explain how it was derived.
This identity is evaluated first by calculating $\frac{\partial I}{\partial a}$ which is
$$
\int -2 a e^{-a^2 \csc ^2(\theta )} \csc ^2(\theta ) \, d\theta = \sqrt{\pi } e^{-a^2} \text{erf}(a \cot (\theta )),
$$
by making the change of variable $\cot \theta = y$.
Next, $I$ is recovered by integrating the RHS of the above equation with regard to $a$ using the fact that this integral amounts to calculating the CDF of skew-normal distribution](https://en.wikipedia.org/wiki/Skew_normal_distribution).