Two electron integral evaluation over primitive Gaussians without recurrence relations
As you mentioned, most implementations such as Obara–Saika (OS), McMurchie–Davidson (MD), Head-Gordon–Pople (HGP), and many others use a series of horizontal and vertical recurrence relations to build up angular momentum and transfer it among each of the four centers. These are analytic solutions, up until the need to solve the Boys function, however we could debate as to whether or not they are closed-form expressions by terminology. They are closed-form in the sense that each part of the recursive expansion does terminate properly with $(ss|ss)$, given in OS equation 44.

In general, any recursive function can be rewritten to be non-recursive using loops of sums and or products, though it may involve considerable effort. In that sense, the recurrence relationships are quite powerful and avoid lots of "boilerplate". If you are interested in a non-recursive expression, there is the earlier Taketa–Huzinaga–O-ohata (THO) paper that, while not very descriptive, includes explicit working expressions for 2-electron and 1-electron (overlap, kinetic energy, and nuclear-electron attraction) integrals. The most relevant expression is

\begin{align*}
\text{ERI} &amp;= \frac{2\pi^{2}}{\gamma_{1}\gamma_{2}} \left( \frac{\pi}{\gamma_{1} + \gamma_{2}} \right)^{1/2} \exp \left\{ - \frac{\alpha_{1}\alpha_{2}\overline{AB}^{2}}{\gamma_{1}} - \frac{\alpha_{3}\alpha_{4}\overline{CD}^{2}}{\gamma_{2}} \right\} \\
&amp;\times \sum_{i_{1}i_{2}r_{1}r_{2}u} B_{i_{1}i_{2}r_{1}r_{2}u}(l_{1}, l_{2}, A_{x}, B_{x}, P_{x}, \gamma_{1} | l_{3}, l_{4}, C_{x}, D_{x}, Q_{x}, \gamma_{2}) \\
&amp;\times \sum_{j_{1}j_{2}s_{1}s_{2}v} B_{j_{1}j_{2}s_{1}s_{2}v}(m_{1}, m_{2}, A_{y}, B_{y}, P_{y}, \gamma_{1} | m_{3}, m_{4}, C_{y}, D_{y}, Q_{y}, \gamma_{2}) \\
&amp;\times \sum_{k_{1}k_{2}t_{1}t_{2}w} B_{k_{1}k_{2}t_{1}t_{2}w}(n_{1}, n_{2}, A_{z}, B_{z}, P_{z}, \gamma_{1} | n_{3}, n_{4}, C_{z}, D_{z}, Q_{z}, \gamma_{2}) F_{\nu}(\overline{PQ}^{2}/4\delta),
\label{tho-2.22-1}\tag{THO eq. 2.22 part 1}
\end{align*}

with

\begin{align*}
&amp;B_{i_{1}i_{2}r_{1}r_{2}u}(l_{1}, l_{2}, A_{x}, B_{x}, P_{x}, \gamma_{1} | l_{3}, l_{4}, C_{x}, D_{x}, Q_{x}, \gamma_{2}) \\
&amp;\quad = (-)^{i_{2}} f_{i_{1}}(l_{1}, l_{2}, \overline{PA_{x}}, \overline{PB_{x}}) f_{i_{2}}(l_{3}, l_{4}, \overline{QC_{x}}, \overline{QD_{x}}) \\
&amp;\qquad \times \frac{i_{1}!i_{2}!}{(4\gamma_{1})^{i_{1}}(4\gamma_{2})^{i_{2}}\delta^{i_{1} + i_{2}}} \cdot \frac{4\gamma_{1})^{r_{1}}(4\gamma_{2})^{r_{2}}\delta^{r_{1} + r_{2}}}{r_{1}!r_{2}!(i_{1}-2r_{1})!(i_{2}-2r_{2})!} \\
&amp;\qquad \times [i_{1} + i_{2} - 2(r_{1} + r_{2})]! \frac{(-)^{u} p_{x}^{i_{1} + i_{2} - 2(r_{1} + r_{2}) - 2u} \delta^{u}}{u! \{ i_{1} + i_{2} - 2(r_{1} + r_{2} - 2u)\}!},
\label{tho-2.22-2}\tag{THO eq. 2.22 part 2}
\end{align*}

and intermediate quantities are

\begin{equation}
\nu = i_{1} + i_{2} + j_{1} + j_{2} + k_{1} + k_{2} - 2 (r_{1} + r_{2} + s_{1} + s_{2} + t_{1} + t_{2}) - u - v - w,
\end{equation}

\begin{align*}
\gamma_{1} &amp;= \alpha_{1} + \alpha_{2},~\gamma_{2}\text{ is similar}\\
\delta &amp;= \frac{1}{4\gamma_{1}} + \frac{1}{4\gamma_{2}} \\
\mathbf{P} &amp;= \frac{\alpha_{1} \mathbf{A} + \alpha_{2} \mathbf{B}}{\alpha_{1} + \alpha_{2}},~\mathbf{Q}\text{ is similar} \\
\mathbf{p} &amp;= \mathbf{Q} - \mathbf{P} \\
\overline{AB} &amp;= A - B
\end{align*}

\begin{equation}
F_{\nu}(t) = \int_{0}^{1} u^{2\nu} \exp(-tu^{2}) \mathop{du},~\text{the Boys function.} \tag{THO eq. 2.11}
\end{equation}

$(l, m, n)$ are the usual polynomial exponents leading to s-, p-, and d-type functions, $(\mathbf{A}, \mathbf{B}, \mathbf{C}, \mathbf{D})$ are the usual basis function positions, the $\alpha$ are the basis function exponents, and $f_{j}(l, m, a, b)$ is the coefficient of $x^{j}$ in the expansion of $(x + a)^{l} (x + b)^{m}$.


  In \eqref{tho-2.22-1} the summations with respect to the indices $i_{1}, i_{2}, r_{1}, r_{2}, u$ extend from 0 to $l_{1} + l_{2}$, $l_{3} + l_{4}$, $[i_{1}/2]$, $[i_{2}/2]$ and $[(i_{1} + i_{2})/2 - r_{1} - r_{2}]$ respectively. The ranges of $(j_{1}, j_{2}, s_{1}, s_{2}, v)$ or $(k_{1}, k_{2}, t_{1}, t_{2}, w)$ may be found in the same way.


Regarding the indices used in the summations:


If the index $i_{1}$ runs from 0 to $l_{1} + l_{2}$, and $i_{2}$ runs from 0 to $l_{3} + l_{4}$ we can deduce that $i_{1}$ must cover all possible values of the total angular momentum along the x-component for the bra, and $i_{2}$ goes up to the total angular momentum along the x-component for the ket. The indices $j,k$ match with $m,n$, meaning the y- and z-component.
$r,s,t$ are not to be confused with distances themselves, like how we usually define $r_{A} = |\mathbf{r}_{A}| = |\mathbf{r} - \mathbf{A}|$. Both $r,s,t$ and $u,v,w$ are also related to the three Cartesian coordinates, but in a more complicated manner than $i,j,k$; they appear to originate from the closed-form expressions for the Hermite polynomials.


Their implementation is slightly modified, with the necessary expressions in equations 3.3, 3.4, and 3.5, though the above can be implemented directly as well. A quick search revealed this code that implements all the THO integrals in C++. It is worth commenting that in real packages, the Boys function is precomputed from interpolation tables, rather than from the gamma and incomplete gamma functions.