Ok, it took me a few thoughtful minutes to unpack @metamorphy's concise answer, but I've finally managed to set things in order and I've decided to collect what I got out of it in a separate answer. @Szeto's formula
$$c_n = \frac{d^n}{dx^n} \frac{2\Gamma(x+1)}{\Gamma(x/2 + 1)^2}\Bigg|_{x=0}$$
can be thought of as providing the $n$-th coefficient of a Maclaurin series, which converges in an open symmetric neighborhood of its center:
$$ \frac{\Gamma(x+1)}{\Gamma(x/2 + 1)^2} = \frac 1 2\sum_{n=0}^\infty \frac{c_n}{n!} x^n .$$
Taking logarithms of both sides (which is allowed because they define nonnegative real numbers on said neighborhood), we obtain
$$\ln\Gamma(x+1)-2\ln\Gamma(x/2+1) = \ln\left[\frac 1 2\sum_{n=0}^\infty \frac{c_n}{n!} x^n\right]. \tag 1$$
Letting $\gamma$ be the Euler-Mascheroni constant, by means of the following log-Gamma series
$$\ln \Gamma(x+1) = -\gamma x + \sum_{k=2}^\infty \frac{(-1)^k \zeta(k)}{k}x^k $$
we may rewrite the LHS of $(1)$ as
$$\ln\Gamma(x+1) - 2\ln\Gamma(x/2 + 1) = \sum_{k=2}^\infty \left[\frac{(-1)^k \zeta(k)}{k} - 2 \frac{(-1)^k \zeta(k)}{k2^k} \right] x^k = \sum_{k=1}^\infty \frac{a_k}{k!} x^k, $$
where we have defined the coefficients
$$a_k := \begin{cases}
(-1)^k (k-1)! (1-2^{1-k}) \zeta(k) = (-1)^k \Gamma(k) \eta(k), & k \geq 2, \\
0, & k = 1,
\end{cases} $$
Going back to $(1)$ and exponentiating both sides, we find
$$\exp\left[ \sum_{k=1}^\infty \frac{a_k}{k!} x^k \right] = \frac 1 2\sum_{n=0}^\infty \frac{c_n}{n!} x^n . \tag2$$
Digression: Bell polynomials. Now the exponential of a power series may be evaluated by resorting to complete Bell polynomials:
$$B_n(x_1,\dots,x_n) = \sum_{k=1}^{n} B_{n,k}(x_1,\dots,x_{n-k+1}) $$
where $B_{n,k}$ is the partial Bell polynomial of order $(n,k)$, defined in the following combinatorial fashion. Set $\mathfrak I_n := \{1,\dots,n\}$ and call $\Pi(n)$ the set of all partitions of the set $\mathfrak I_n$; so that,
$$\begin{split}
p \in \Pi(n) \quad \iff \quad &p = \{S_1, \dots S_k\},\ \text{for some } k \in \mathfrak I_n,\ \text{s.t.} \\ &\forall \ell \ \ \varnothing \neq S_\ell \subseteq \mathfrak I_n; \quad \bigcup_{\ell=1}^k S_\ell = \mathfrak I_n; \quad S_\ell \cap S_m = \varnothing.
\end{split}$$
Now let $k \in \mathfrak I_n$. Then there is a certain number $s_{n,k}$ of partitions of $\mathfrak I_n$ into $k$ blocks (the Stirling number of the second kind, also denoted ${n\brace k}$). These blocks can be of various sizes (at most of size $n-k+1$ as empty blocks are forbidden), obviously adding up to $n$. But only a certain amount of those partitions will be such that the size $ \sigma(\ell)$ of the $\ell$-th block, for all $\ell \in \mathfrak I_k$, is a fixed number. Let $\sigma : \mathfrak I_k \to \mathfrak I_{n-k+1}$ denote such an assignment of sizes, and let the amount of partitions of $\mathfrak I_n$ into $k$ blocks with the mandated sizes be $b_{n,k,\sigma}$; if $\Sigma(n,k)$ is the collection of all possible $\sigma$'s (taken modulo permutations of the same sizes among the blocks), then $B_{n,k}$ is the polynomial
$$B_{n,k}(x_1,\dots,x_{n-k+1}) = \sum_{\sigma \in \Sigma(n,k)} b_{n,k,\sigma} \prod_{\ell = 1}^k x_{\sigma(\ell)}$$
Since as we said the total number of partitions into $k$ parts must be ${n\brace k}$, we have for instance
$$B_{n,k}(1,\dots,1) = \sum_{\sigma \in \Sigma(n,k)} b_{n,k,\sigma} = {n\brace k}; $$
summing over all possible values of $k$, we obtain
$$B_n(1,\dots,1) = \sum_{k=1}^n B_{n,k}(1,\dots,1) = \sum_{k=1}^n {n\brace k} =: B_n, $$
which defines the $n$-th Bell number. As for the coefficients $b_{n,k,\sigma}$, they may be calculated in the following way: first, suppose $\sigma$ is injective, so that each block in the partitions it defines has a different size than the others. Then
$$b_{n,k,\sigma} = \binom{n}{\sigma(1)} \binom{n-\sigma(1)}{\sigma(2)} \cdots \binom{n-\sigma(1) - \dots - \sigma(k-2)}{\sigma(k-1)} \binom{n-\sigma(1) - \dots - \sigma(k-1)}{\sigma(k)}, $$
because there are $\binom{n}{\sigma(1)}$ ways to choose $\sigma(1)$ objects out of a collection of $n$ objects, and $\binom{n-\sigma(1)}{\sigma(2)}$ ways to choose $\sigma(2)$ objects out of the remaining $n-\sigma(1)$ objects, and so on. This simplifies to
$$b_{n,k,\sigma} = \frac{n!}{\sigma(1)! \cdots \sigma(k)!} = \binom{n}{\sigma(1),\dots,\sigma(k)}. $$
If $\sigma$ fails to be injective, so that some blocks have the same size as others, the formula above ends up overcounting the number of possibilities (because it treats blocks with the same size on different grounds). Then it must be divided out, for each repeated size, by the factorial of the total number of blocks that have that size – which is the total number of ways to permute those blocks between each other without changing the setup of the partition. To do so, we must have knowledge of where and how badly the map $\sigma$ fails to be injective: this can be done by looking at the matrix
$$[\sigma]=\begin{bmatrix}
\sigma(1) & & \\
& \ddots & \\
& & \sigma(k)
\end{bmatrix} \in \operatorname{Mat}_k(\mathbb Z_{n-k+1}) $$
and analyzing its algebraic degeneracy.
A polynomial form for $c_n$. Going back to $(2)$, one may show that
$$\exp\left[ \sum_{k=1}^\infty \frac{a_k}{k!} x^k \right] = \sum_{n=0}^\infty \frac{B_n(a_1,\dots,a_n)}{n!}x^n, $$
so that we may compare RHS's:
$$\frac 1 2\sum_{n=0}^\infty \frac{c_n}{n!} x^n = \sum_{n=0}^\infty \frac{B_n(a_1,\dots,a_n)}{n!}x^n, \quad \implies \quad \boxed{ c_n = 2B_n(a_1,\dots,a_n).}$$
