How many expressions can be formed with $k$ commutative and associative functions? 
This is a generalization of a question I posted to Math.SE here.

Suppose we have $k$ binary functions $f_1, f_2, \ldots, f_k$, all of which are commutative and associative, and $n$ indeterminates $x_1, x_2, \ldots, x_n$. What is the number $a_{n,k}$ of distinct expressions can we form using any combination (repetitions allowed) of the $f_i$ and one each of the $x_i$? Is there any efficient way to generate all of these expressions?
For example, in the case $k=2$, $n=3$, we can form $a_{3,2} = 8$ distinct expressions:
$$f_1(f_1(x_1,x_2),x_3) \quad f_1(f_2(x_1,x_2),x_3) \quad f_1(f_2(x_1,x_3),x_2) \quad f_1(f_2(x_2,x_3),x_1)$$
$$f_2(f_1(x_1,x_2),x_3) \quad f_2(f_1(x_1,x_3),x_2) \quad f_2(f_1(x_2,x_3),x_1) \quad f_2(f_2(x_1,x_2),x_3)$$
After poking around on the OEIS for a while, I conjecture (for reasons entirely beyond my understanding) that
$$ a_{n,k} = \sum_{i=0}^n \sum_{j=0}^i (-1)^j k^i {i+n-1 \choose j+n-1} {j+n-1 \brack j} $$
where ${n \brack k}$ is a Stirling number of the first kind. I have checked this for $n,k = 1$ to $4$, but I have no idea if or why it works.
 A: One way to formulate the (first) question you have asked is:

What is the dimension of the arity $n$ component of the free product of $k$ copies of the commutative operad?

Koszul duality for operads gives the following method of calculation (spoiler: I'm not going to end with an explicit formula). Fix $k$ and then combine the terms in an exponential generating function:
$$
Q_k(x) = \sum_{n=1}^\infty \frac{a_{n,k}}{n!}x^n.
$$
In order to calculate $Q_k(x)$ we'll examine another, seemingly unrelated question which is easier to answer: Suppose we have $k$ binary functions which are mutually annihilating Lie brackets; then how many distinct expressions can we make from $n$ letters? For $k=1$ the answer is classically known to be $(n-1)!$; then since they are mutually annihilating, the answer in the general case is
$$
b_{n,k} = \left\{
\begin{array}{ll}
1&n=1\\
k(n-1)!&\text{otherwise.}
\end{array}
\right.
$$
Combining these terms into a exponential generating series as above we get
$$
P_k(x) = \sum_{n=1}^\infty \frac{b_{n,k}}{n!}x^n = x+\sum_{n=2}^\infty \frac{k}{n}x^n = (1-k)x -k\log(1-x).
$$
Here we have calculated the exponential generating function for the arity $n$ component of the direct sum of $k$ copies of the Lie operad.
Then classical results about Koszul duality theory for operads assert the following:


*

*the commutative and Lie operads are Koszul dual,

*if $A$ and $B$ are Koszul dual to $A'$ and $B'$, then the free product of $A$ and $B$ is Koszul dual to the direct sum of $A'$ and $B'$, and

*if two operads are Koszul dual, then their exponential generating series $P(x)$ and $Q(x)$ are related by $-P(-x)=Q^{-1}(x)$ (in formal series expansions near zero).


Together this tells us that $Q_k(x)$ is the formal inverse of the function
$$
Q_k^{-1}(x) = (1-k)x +k\log(1+x).
$$
Of course, $Q_k(x)$ is not elementary but it can be written in terms of the Lambert $W$ function. It's something like the following (but I think I might have made a typo or a sign error, check this formula!)
$$
Q_k(x)=-1-\frac{k}{k-1} W\left(\frac{k-1}{k}e^{\frac{x+1-k}{k}}\right).
$$
In any event, one can use purely symbolic analysis of the $W$ function and a single special evaluation 
$$W\left(\frac{1-k}{k}e^{\frac{1-k}{k}}\right)=\frac{1-k}{k}
$$
 to work out a series expansion. 
For instance, for $k=2$ we get:
$$
x + x^2 + \frac{4x^3}{3} + \frac{13x^4}{6} + \frac{59x^5}{15}+\cdots
$$
which gives
$a_{1,2}=1$
$a_{2,2}=2$
$a_{3,2}=8$
$a_{4,2}=52$
$a_{5,2}=472$.
I'm too lazy to check whether this agrees with anything in the OEIS but this series (and those for other $k$), presumably, is well-known and can be cross-referenced there.
