How many directed graphs of size n are there where each vertex is the tail of exactly one edge?

Here are some computational aspects of this problem that can serve as a basis for further exploration. We present relevant generating functions and recurrence relations for unlabelled trees and endofunctions. These basic data of course cannot take the place of a proper investigation using combinatorial methods.

We have by inspection that the combinatorial class $\mathcal{M}$ under consideration here is (notation from Analytic Combinatorics)

$$\mathcal{M} = \def\textsc#1{\dosc#1\csod} \def\dosc#1#2\csod{{\rm #1{\small #2}}}\textsc{MSET}(\textsc{CYC}(\mathcal{T}))$$ with $\mathcal{T}$ being the combinatorial class of rooted trees. Therefore the study of $\mathcal{M}$ requires the study of $\mathcal{T}.$

We have for rooted trees that $$\mathcal{T} = \mathcal{Z} \times \textsc{MSET}(\mathcal{T}).$$ This translates to the functional equation $$T(z) = z \exp\left(\sum_{l\ge 1} \frac{T(z^l)}{l} \right).$$

Differentiate to get a recurrence for the number $T_n$ of rooted trees: $$T'(z) = \exp\left(\sum_{l\ge 1} \frac{T(z^l)}{l} \right) + z \exp\left(\sum_{l\ge 1} \frac{T(z^l)}{l} \right) \left(\sum_{l\ge 1} \frac{T'(z^l)}{l} \times l z^{l-1} \right) \\ = \frac{T(z)}{z} + T(z) \left(\sum_{l\ge 1} T'(z^l) z^{l-1} \right) .$$ so that $$z T'(z) = T(z) + T(z) \left(\sum_{l\ge 1} T'(z^l) z^l \right) = T(z) + T(z) \left(\sum_{l\ge 1} z^l \sum_{m\ge 1} m T_m z^{(m-1)l}\right)$$ which finally yields $$z T'(z) = T(z) + T(z) \left(\sum_{l\ge 1} \sum_{m\ge 1} m T_m z^{ml}\right).$$

Extracting coefficients from this we get $$n T_n = T_n + [z^n] T(z) \left(\sum_{l\ge 1} \sum_{m\ge 1} m T_m z^{ml}\right) = T_n + \sum_{l=1}^{n-1} \sum_{m=1}^{\lfloor(n-1)/l\rfloor} m T_m T_{n-ml}$$ so that $$T_n = \frac{1}{n-1} \sum_{l=1}^{n-1} \sum_{m=1}^{\lfloor(n-1)/l\rfloor} m T_m T_{n-ml}$$ with $T_0 = 0$ and $T_1 = 1.$

This gives the sequence $$1, 1, 2, 4, 9, 20, 48, 115, 286, 719, 1842, 4766, 12486, 32973, 87811, \ldots$$ which is OEIS A000081 where the above derivation is confirmed.

Now return to the combinatorial class $\mathcal{Q}$ and compute the generating function of the operator $\textsc{MSET}(\textsc{CYC}(\cdot)).$

Using the cycle index of the cyclic group we easily derive that the generating function of the cycle operator being applied to some combinatorial class $\mathcal{A}$ is $$\sum_{k\ge 1} \frac{\varphi(k)}{k} \log\frac{1}{1-A(z^k)}.$$

Hence the generating function for $\textsc{MSET}(\textsc{CYC}(\cdot))$ is $$\exp\left( \sum_{l\ge 1}\frac{1}{l} \sum_{k\ge 1} \frac{\varphi(k)}{k} \log\frac{1}{1-A(z^{kl})}\right).$$ This is $$\exp \left(\sum_{m\ge 1} \log\frac{1}{1-A(z^m)} \sum_{k|m} \frac{\varphi(k)}{m}\right) \\ = \exp \left(\sum_{m\ge 1} \log\frac{1}{1-A(z^m)}\right) = \prod_{m\ge 1} \frac{1}{1-A(z^m)}.$$

It follows that the generating function $Q(z)$ of $\mathcal{Q}$ in terms of $T(z)$ is given by $$Q(z) = \prod_{m\ge 1} \frac{1}{1-T(z^m)}.$$

Differentiate to obtain a recurrence from this: $$Q'(z) = \prod_{m\ge 1} \frac{1}{1-T(z^m)} \times \sum_{m\ge 1} (1-T(z^m)) \frac{T'(z^m) \times m z^{m-1}}{(1-T(z^m))^2} \\ = Q(z) \times \sum_{m\ge 1} \frac{T'(z^m) \times m z^{m-1}}{1-T(z^m)}.$$

Extracting coefficients from this we obtain $$(n+1) Q_{n+1} = \sum_{q=0}^n Q_{n-q} [z^q] \sum_{m\ge 1} \frac{T'(z^m) \times m z^{m-1}}{1-T(z^m)} \\ = \sum_{q=0}^n Q_{n-q} [z^q] \sum_{m\ge 1} \frac{m}{z} \frac{\sum_{p\ge 1} p T_p z^{m (p-1)} \times z^m}{1-T(z^m)} \\ = \sum_{q=0}^n Q_{n-q} \sum_{m=1}^{q+1} m [z^{q+1}] \frac{\sum_{p\ge 1} p T_p z^{mp}}{1-T(z^m)}.$$

Now here we need to digress for a moment to find a recurrence for the coefficients $P_n$ of $$P(z) = \frac{1}{1-T(z)}.$$ Differentiating yields $$P'(z) = \frac{T'(z)}{(1-T(z))^2} = P(z)^2 T'(z).$$ Extracting coefficients we obtain $$(n+1) P_{n+1} = \sum_{q=0}^n \left(\sum_{p=0}^{n-q} P_p P_{n-q-p}\right) [z^q] T'(z) \\ = \sum_{q=0}^n (q+1) T_{q+1} \left(\sum_{p=0}^{n-q} P_p P_{n-q-p}\right)$$ so that $$P_n = \frac{1}{n} \sum_{q=0}^{n-1} (q+1) T_{q+1} \left(\sum_{p=0}^{n-1-q} P_p P_{n-1-q-p}\right)$$ with $P_0 = 1.$

This yields the sequence $$1, 2, 5, 13, 35, 95, 262, 727, 2033, 5714, 16136, 45733, 130046, 370803,\ldots$$ which is OEIS A000107 where a better recurrence can also be found.

Returning to $Q_{n+1}$ and using an Iverson bracket we obtain $$\sum_{q=0}^n Q_{n-q} \sum_{m=1}^{q+1} m [[m|q+1]] \sum_{p=1}^{(q+1)/m} p T_p P_{(q+1)/m-p}.$$ which yields the formula $$Q_n = \frac{1}{n} \sum_{q=0}^{n-1} Q_{n-1-q} \sum_{m=1}^{q+1} m [[m|q+1]] \sum_{p=1}^{(q+1)/m} p T_p P_{(q+1)/m-p}$$ where $Q_0 = Q_1 = 1.$

An alternate form of this is $$Q_n = \frac{1}{n} \sum_{q=0}^{n-1} Q_{n-1-q} \sum_{m|q+1} m \sum_{p=1}^{(q+1)/m} p T_p P_{(q+1)/m-p}$$ again where $Q_0 = Q_1 = 1.$

This finally gives the sequence $$1, 3, 7, 19, 47, 130, 343, 951, 2615, 7318, 20491, 57903, 163898, \\ 466199, 1328993,\ldots$$ which points us to OEIS A001372 where additional material awaits and the results of the above calculation are validated.

The following Maple code implements these three sequences.

with(numtheory, divisors);

T :=
proc(n)
option remember;

    if n<=1 then return n fi;

    1/(n-1)*add(add(m*T(m)*T(n-m*l),
                    m=1..floor((n-1)/l)), l=1..n-1);
end;

P :=
proc(n)
    option remember;

    if n=0 then return 1 fi;

    1/n*add((q+1)*T(q+1)*
            add(P(p)*P(n-1-q-p), p=0..n-1-q), q=0..n-1);

end;

Q :=
proc(n)
option remember;
    local res, iv, q, m;

    if n<=1 then return 1 fi;

    1/n*add(Q(n-1-q)*
            add(m*add(p*T(p)*P((q+1)/m-p), p=1..(q+1)/m),
                m in divisors(q+1)), q=0..n-1);
end;

Addendum. Note that just for curiosity's sake we can do the labelled case es well, which is very simple. Labelled trees are given by $$\mathcal{T} = \mathcal{Z} \times \textsc{SET}(\mathcal{T})$$ which gives the functional equation $$T(z) = z \exp T(z).$$

The combinatorial class of endofunctions is then given by $$\mathcal{Q} = \textsc{SET}(\textsc{CYC}(\mathcal{T})).$$

Translating to generating functions we get $$Q(z) = \exp \log \frac{1}{1-T(z)} = \frac{1}{1-T(z)}.$$

Extracting coefficients via Lagrange inversion we have $$Q_n = n! \frac{1}{2\pi i} \int_{|z|=\epsilon} \frac{1}{z^{n+1}} \frac{1}{1-T(z)} dz.$$

Put $T(z)=w$ so that $z=w/\exp(w) = w\exp(-w)$ and $dz = \exp(-w) - w\exp(-w)$ to get $$n! \frac{1}{2\pi i} \int_{|w|=\epsilon} \frac{\exp(w(n+1))}{w^{n+1}} \frac{1}{1-w} (\exp(-w) - w\exp(-w)) dw \\ = n! \frac{1}{2\pi i} \int_{|w|=\epsilon} \frac{\exp(wn)}{w^{n+1}} \frac{1}{1-w} (1 - w) dw \\ = n! \frac{1}{2\pi i} \int_{|w|=\epsilon} \frac{\exp(wn)}{w^{n+1}} dw.$$

But we have $$n! [w^n] \exp(wn) = n! \times \frac{n^n}{n!} = n^n$$

and we have just verified the count of the number of endofunctions from basic algebra.

All of the above (labelled and unlabelled) is already present in some form or other in Harary and Palmer's Graphical Enumeration.

Remark I. Care should be taken not to confuse the labelled and unlabelled tree functions.

Remark II. Care should furthermore be taken in the proper usage of sets and multisets in the labelled and the unlabelled case. We distinguish multisets and sets in the unlabelled case but in the labeled case the cartesian product that underlies the construction of all combinatorial operators involves re-labeling that uniquely identifies the two components by the choice of labels and there can never be multisets (operator $\textsc{MSET}$), just sets (operator $\textsc{SET}$).

A beginning:

Let $f:\>V\to V$ be the map defined by $f(v):=$ endpoint of the edge emanating from $v\in V$.

For an intuitive understanding of what's going on here you should consider the iterates $f^{\circ k}$, $\>k\geq0$. Since $V$ is finite each orbit $$O(v):=\{f^{\circ k}(v)\>|\>k\geq0\}$$ must end in a cycle. Call two vertices equivalent if their orbits end in the same cycle. The equivalence classes, together with the connecting edges, are the components of your graph $\Gamma$.

Such a component has the following shape: There is a final cycle $\gamma=\{v_1,v_2,\ldots, v_r\}$ of length $r\geq1$ if loops are allowed in $\Gamma$, and $r\geq2$ otherwise. Each $v_i\in\gamma$ is the root of a rooted tree collecting the vertices whose orbits enter $\gamma$ at $v_i$.

I hope that such things can be counted up to isomorphism using Polya counting theory. (The counting of rooted trees is a standard example in this theory. But here the situation is more complicated, since we have necklaces of rooted trees.)