I was working on this question and come up with this problem.
Let $\sigma_i,\sigma_j$ be two uniformly chosen derangements, i.e. $\sigma_i,\sigma_j \in D_n = \{\sigma \in S_n : \sigma(i)\neq i~\forall i\}$. What is the distribution of $\sigma_i\circ \sigma_j$?
Here we can see that every $\sigma \in S_n$ is the composition of two $\sigma_i,\sigma_j\in D_n$. But some permutations are easier to produce than others. For instance, let $\sigma^*$ be the identity, we have $\sigma_i\circ\sigma_j = \sigma^*$ iif $\sigma_j =\sigma_i^{-1}$. So, we are free to choose $\sigma_i$, but then $\sigma_j$ is defined, therefore: $$P(\sigma_i\circ\sigma_j = \sigma^*) = \frac{1}{|D_n|} > \frac{1}{|S_n|} $$
But $\lim_{n\to\infty}\frac{|D_n|}{|S_n|} = \frac{1}{e}$ (see here), it leads me to think that although it is not uniform, it could at least have the same order. So, my guess is that we can find $c,C > 0$ such that $$c\frac{1}{|S_n|} \leq P(\sigma_i\circ\sigma_j = \sigma) \leq C \frac{1}{|S_n|}~\forall \sigma\in S_n$$ for sufficiently large values of $n$.
I'll find an approximation for the distribution below, which proves your guess is right. Note that for any given $\pi \in S_n$, $\mathbb{P}[\sigma_i \circ \sigma_j = \pi]$ is $\frac{1}{|D_n|^2}$ times the number $N(\pi)$ of derangements $\sigma$ such that $\pi^{-1} \circ \sigma$ is also a derangement. Equivalently, $N(\pi)$ is the number of $\sigma \in S_n$ such that $\sigma(x) \neq x, \pi(x)$ for all $x \in [n]$. We can find an expression for $N(\pi)$ using an inclusion-exlusion argument similar to the standard one used to find the number of derangements. We have $$N(\pi) = \sum_{A \subset [n]} (-1)^{|A|} \#\{\sigma \in S_n \mid \sigma(a) \in \{a, \pi(a)\} \, \forall a \in A\}.$$ Now, given a set $A \subset [n]$ and an injective function $f : A \to [n]$ such that $f(a) \in \{a, \pi(a)\}$ for all $a \in A$, there are exactly $(n-|A|)!$ choices of $\sigma \in S_n$ which agree with $f$ on $A$, so we can write $$\#\{\sigma \in S_n \mid \sigma(a) \in \{a, \pi(a)\} \, \forall a \in A\} = F_A(\pi)(n-|A|)!$$ where $F_A(\pi)$ is the number of such functions. This gives $$N(\pi) = \sum_{A \subset [n]} (-1)^{|A|} F_A(\pi)(n-|A|)! = n! \sum_{k=0}^n \frac{(-1)^k}{k!} \binom{n}{k}^{-1} \sum_{|A| = k} F_A(\pi).$$ Clearly each $F_A(\pi) \leq 2^k$, so $\binom{n}{k}^{-1} \sum_{|A| = k} F_A(\pi) \leq 2^k$, so these terms decrease extremely quickly: $$\left|\frac{N(\pi)}{n!} - \sum_{k=0}^m \frac{(-1)^k}{k!} \binom{n}{k}^{-1} \sum_{|A| = k} F_A(\pi) \right| \leq \sum_{k=m+1}^n \frac{2^k}{k!} \leq \frac{2^m}{m!} \sum_{k=1}^\infty \frac{2^k}{m^k} \leq \frac{2^m}{m!}$$ holds for $m \geq 4$.
Now we find an approximation for the largest terms. Note that $k! \sum_{|A| = k} F_A(\pi)$ is the number of sequences of pairs $(a_1, b_1), \dots, (a_k, b_k) \in [n]^2$ such that all $a_i$ are distinct, all $b_i$ are distinct, and $b_i \in \{a_i, \pi(a_i)\}$ for each $i$. Call a sequence $a_1, \dots, a_k$ good if all $\{a_i, \pi(a_i)\}$ are disjoint (so in particular all $a_i$ are distinct), and note that for a good sequence, the number of corresponding sequences of pairs is $\prod_{i=1}^k |\{a_i, \pi(a_i)\}|$. As before, for any (not necessarily good) sequence $a_1, \dots, a_k$, the number of corresponding sequences of pairs is at most $2^k$. Finally, the number of sequences $a_1, \dots, a_k$ which are not good is at most $2k^2 n^{k-1}$ (since a sequence is not good if there are indices $i$ and $j$ with either $a_i = a_j$ or $a_i = \pi(a_j)$). Thus $$\left|k!\sum_{|A| = k}F_A(\pi) - \sum_{a_1, \dots, a_k} \prod_{i=1}^k |\{a_i, \pi(a_i)\}|\right| \leq 2^{k+1}k^2 n^{k-1}$$ but $\sum_{a_1, \dots, a_k} \prod_{i=1}^k |\{a_i, \pi(a_i)\}| = (\sum_a |\{a, \pi(a)\}|)^k = (2n-t)^k$, where $t$ is the number of fixed points of $\pi$. Letting $c = t/n$, we have $\left|\frac{k!}{n^k} \sum_{|A| = k}F_A(\pi) - (2-c)^k \right| \leq \frac{2^{k+1}k^2}{n}.$ This can be written as $\frac{k!}{n^k} \sum_{|A| = k}F_A(\pi) = (2-c)^k + O(\frac{2^k k^2}{n})$, and so since $\binom{n}{k}^{-1} = \frac{k!}{n^k}(1 + O(\frac{k^2}{n}))$ for small $k$, we get $$\binom{n}{k}^{-1} \sum_{|A| = k}F_A(\pi) = (2-c)^k + O\left(\frac{2^k k^2}{n}\right)$$ hence $\sum_{k=0}^m (-1)^k \binom{n}{k}^{-1} \sum_{|A| = k}F_A(\pi) = \sum_{k=0}^m \frac{(c-2)^k}{k!} + O(1/n)$. Noting that $|\sum_{k=0}^\infty \frac{(c-2)^k}{k!} - \sum_{k=0}^m \frac{(c-2)^k}{k!}| \leq \frac{2^m}{m!}$ for $m \geq 4$ as well, and taking $m$ large enough so that $2^m/m! \leq 1/n$, we get $$\frac{N(\pi)}{n!} = \sum_{k=0}^\infty \frac{(c-2)^k}{k!} + O(1/n) = e^{c-2} + O(1/n)$$ where the error term is uniform over all $\pi$. Thus since $|D_n| = n! e^{-1}(1 + O(1/n))$, we have $$\mathbb{P}[\sigma_i \circ \sigma_j = \pi] = \frac{N(\pi)}{|D_n|^2} = \frac{1}{n!} (e^c + O(1/n))$$ where again, $c = c(\pi) = t/n$ is the fraction of $[n]$ fixed by $\pi$.