Closed form of $\int_0^1 \tan(\gamma\sqrt{1-x^2}) dx$

Some context: I'm studying the problem of nonperturbative pair creation from strong fields in quantum electrodynamics. For certain time dependent electric fields I can get some information about the qualitative behavior from integrals of the form $$ \tilde{h}(\gamma) = \int_0^1 h(\gamma\sqrt{1-x^2}) dx, $$ where $h$ depends on the field profile and $\gamma$ is something like an inverse timescale (if $h$ has any poles, $\gamma$ is always restricted to an appropriate region).

For $h=\sinh$, $h=\tanh^{-1}$ and $h=\mathrm{erfi}$ (the imaginary error function) Mathematica provides closed form solutions, but not for $h=\tan$, which I am greatly interested in.

The only idea I had was expressing $h$ using its Taylor series and integrating term by term, which I can do because

$$ \int_0^1 \gamma^n (1-x^2)^{n/2} dx = \gamma^n \frac{\sqrt{\pi}}{2} \frac{\Gamma(1+n/2)}{\Gamma(3/2 + n/2)}, $$

and then sum up the resulting series. This again works for the other functions mentioned above, but not for the tangent.

Any ideas how to get a closed form result of the integral for $h=\tan$?

Edit: The results Mathematica is able to compute are: $$ \text{for } h = \sinh: \quad\tilde{h}(\gamma) = \frac{\pi}{2} I_1(\gamma ), $$ $$ \text{for } h = \tanh^{-1}: \quad\tilde{h}(\gamma) = \frac{\pi}{2}\frac{1-\sqrt{1-\gamma ^2}}{ \gamma }, $$ $$ \text{for } h = \mathrm{erfi}: \quad\tilde{h}(\gamma) = \frac{\sqrt{\pi }\gamma}{2} e^{\frac{\gamma ^2}{2}} \left(I_0\left(\frac{\gamma ^2}{2}\right)-I_1\left(\frac{\gamma ^2}{2}\right)\right). $$

I don't think you will find a nice closed form in the case of $I(\gamma)=\int_0^1 \tan(\gamma\sqrt{1-x^2}) dx$.

I have managed to integrate the series expansion for $\tan(\gamma\sqrt{1-x^2})$ term by term using mathematica. Due to the limits used a lot of the resultant terms are zero. The only term that is non-zero is a term involving $\arcsin(x)$. The result I find is similar in form to @nemo's second result, although nemo uses the falling factorial in the form $\frac{(1/2)_n}{n!}$ which differs from my formulation as it alternates in sign. Assuming I have calculated it correctly, nemo's second result does not agree numerically with the original integral or my result.

My result is

$$I(\gamma)=\frac{\pi}{\gamma}\sum_{k=1}^\infty \left(2^{2k}-1\right)\zeta(2k)\left(\frac{\gamma}{\pi} \right)^{2k}\prod_{j=1}^k\left(\frac{2j-1}{2j}\right),$$

which is different to any of the results given above in the comments.

Substituting for $\zeta(2k)$ using $\zeta(2k)=\frac{(2\pi)^{2k}}{2(2k)!}|B_{2k}|$ and for the product using the identity $\prod_{j=1}^k\left(\frac{2j-1}{2j}\right)=\frac{(2k)!}{2^{2k}(k!)^2} =\frac{1}{2^{2k}}\binom{2k}{k}$ gives

$$I(\gamma)=\frac{\pi}{2\gamma}\sum_{k=1}^\infty \frac{(2^{2k}-1)\;|B_{2k}|\;{\gamma}^{2k}}{(2k)!} \binom{2k}{k}.$$

The rational product $\prod_{j=1}^k\left(\frac{2j-1}{2j}\right)$ could also be written in terms of the ratio of two gamma functions (after using the substitution $\Gamma(1/2)=\sqrt{x}$). This is the one Mathematica likes to use.

$$\prod_{j=1}^k\left(\frac{2j-1}{2j}\right)=\frac{\Gamma \left(k+\frac{1}{2}\right)}{\sqrt{\pi } \;\Gamma (k+1)}.$$

This is interesting as it appears the series for $I(\gamma)$ is a hybrid between two different classes of standard series (e.g. a sort of hybrid between $\tan x$ or $\log \cos x$ involving Bernoulli numbers and $\arcsin x$ involving $\frac{1}{2^{2k}}\binom{2k}{k}$).