By using the exponential generating function for the Legendre polynomials:
\begin{equation}
e^{xz}J_{0}\left(z\sqrt{1-x^{2}}\right)=\sum_{n=0}^{\infty}\frac{P_{n}\left(x\right)}{n!}z^{n}
\end{equation}
with $z=-t,x=\cos\theta$ the first integral (taken with $a=1$, as remarked by @Gary in a comment) reads
\begin{align}
\varphi (R) &= \int_0^\frac{\pi}{2} \int_0^R e^{-t\cos\theta} \sin^3\theta\, J_0(t\sin\theta) \, dt \, d\theta \\
&=\sum_{n=0}^{\infty}\frac{(-1)^n}{n!}\int_0^\frac{\pi}{2} \int_0^Rt^nP_n(\cos\theta)\sin^3\theta \, dt \, d\theta\\
&=\sum_{n=0}^{\infty}\frac{(-1)^n}{(n+1)!}R^{n+1}\int_0^1(1-x^2)P_n(x)\,dx
\end{align}
we supposed valid the interchange of integration and summation and change $x=\cos\theta$ to obtain the latter expression.
When $n$ is even ($n=2p$), we have, by symmetry,
\begin{align}
w_{2p}&=\int_0^1(1-x^2)P_{2p}(x)\,dx\\
&=\frac{1}{2}\int_{-1}^1(1-x^2)P_{2p}(x)\,dx
\end{align}
The orthogonality properties of the Legendre polynomials indicate that the integrals with $p\ge2$ vanish and direct calculation shows that $w_0=2/3,w_2=-2/15$.
When $n$ is odd ($n=2p+1$), we use the tabulated integral (DLMF)
\begin{equation}
\int_{0}^{1}P_{2p+1}\left(x\right)x^{z-1}\,dx=\frac{(-1)^{p}{\left(1-\frac{1}{2}z\right)_{p}}}{2{\left(\frac{1}{2}+\frac{1}{2}z\right)_{p+1}}}
\end{equation}
where $(.)_n$ is Pochhammer's symbol, to express
\begin{align}
w_{2p+1}&=\int_0^1(1-x^2)P_{2p+1}(x)\,dx\\
&=\frac{(-1)^{p}{\left(\frac{1}{2}\right)_{p}}}{2{\left(1\right)_{p+1}}}-\frac{(-1)^{p}{\left(-\frac{1}{2}\right)_{p}}}{2{\left(2\right)_{p+1}}}
\end{align}
Then, converting the series into generalized hypergeometric functions (GHF),
\begin{align}
\varphi(R)&=\sum_{n=0}^{\infty}\frac{(-1)^n}{(n+1)!}R^{n+1}w_n\\
&=w_0R+\frac16w_2R^3-\sum_{p=0}^\infty\frac{1}{(2p+2)!}w_{2p+1}R^{2p+2}\\
&=\frac23R-\frac1{45}R^3-\frac12\sum_{p=0}^\infty\frac{(-1)^p}{(2p+2)!}R^{2p+2}\left[\frac{\left(\frac{1}{2}\right)_{p}}{{\left(1\right)_{p+1}}}-\frac{\left(-\frac{1}{2}\right)_{p}}{{\left(2\right)_{p+1}}}\right]\\
&=\frac23R-\frac1{45}R^3-\frac{R^2}4\left[
\phantom{}_{2}F_{3}\! \left(\frac{1}{2},1;\frac{3}{2},2,2;-\frac{R^{2}}{4}\right)
-\frac12\phantom{}_{2}F_{3}\! \left(-\frac{1}{2},1;\frac{3}{2},2,3;-\frac{R^{2}}{4}\right)
\right]
\end{align}
Then the asymptotic behavior follows:
\begin{equation}
\varphi(R)=\frac23+\sqrt{\frac{2}{\pi}}\frac{\sin(R-\pi/4)}{R^{3/2}}-\frac{21}{8}\sqrt{\frac{2}{\pi}}\frac{\sin(R+\pi/4)}{R^{5/2}}+O(R^{-7/2})
\end{equation}
This result was obtained with a CAS, however it should not be too difficult to derive it from the asymptotic expansion of the GHF.
This method should also apply to evaluate the second integral, by differentiating the starting identity
\begin{equation}
J_0(t\sin\theta)=e^{t\cos(\theta)}\sum_{n\ge0}\frac{(-1)^n}{n!}P_n(\cos\theta)t^n
\end{equation}
w.r.t. $t$ or $\theta$ to obtain an expansion of $J_1$, but I didn't try it.