Constructing Very "Flat" Functions Consider the function $f(x)=\tan(\sin(x))-\sin(\tan(x))$. It has (what I think is) a remarkable property. The Taylor series expansion of $f(x)$ around $x=0$ is
$$\frac{1}{30}x^7+\frac{29}{756}x^9+\frac{1913}{75600}x^{11}+o(x^{13}).$$
Since $\tan(x)$ and $\sin(x)$ are odd functions, $f(x)$ is also an odd function. But its first nonzero term is $x^7$. This tells us that that $f(x)$ is very "flat" near the origin, which it is indeed.

I haven't calculated an explicit closed form for the series coefficients but looking at the wild oscillations around $x=\pm \pi/2$ (and confirming it by looking at it in the complex plane), the radius of convergence is $\pi/2$. This is simply because of the tangent function.
I have the following questions.
Is this just a coincidence or is there something deeper going on here? Is there any way to build/find more such examples? How can we find two real functions $g_1(x)$ and $g_2(x)$ such that

*

*they are both distinct,

*they are both elementary and can be explicitly written in closed form,

*and $g_1(g_2(x))-g_2(g_1(x))$ is $o(x^n)$ for some $n>1$ with a strictly positive radius of convergence.

Neither should be a polynomial because this is trivial to do with polynomials. So the Taylor expansion of $g_1(x)$, $g_2(x)$, and $g_1(g_2(x))-g_2(g_1(x))$ should be infinite.

Addressing Calvin Khor's comment
Is it really trivial with polynomials? The first few terms of the Taylor series only depend on the first few derivatives, so it seems to me that understanding how to do it for polynomials should help understand what's happening in your example
Yes, it is trivial with polynomials. Consider
$$g_1(x)=x^{11}+x^{13}$$
$$g_2(x)=x^{15}+x^{17}$$
then
$$g_2(g_1(x))-g_1(g_2(x))=4 x^{167}+50x^{169}+290x^{171}+\cdots$$
which is indeed a finite series with an infinite radius of convergence, because it is only a composition of polynomials. So it is trivial to build such a (finite series) example for any desired $o(x^n)$.
I was looking for a general construction for any given $n$, so that the series is infinite and can be written in explicit closed form using elementary functions, or a proof that it is not possible.
 A: Summary
$O(x^7)$ is only one power more surprising than most $O(x^3)$ analytic perturbations of the trivial pair $g_1(x)=g_2(x)=x$, which give $O(x^6)$ errors. Some distinct pairs $g_1,g_2$ give exact cancellation. The extra power is explained by the fact that $\sin$ and $\tan$ are odd functions (implying the vanishing of the $x^2$ and $x^4$ terms in their expansion.)
Perfect pairs
There is a pair of distinct elementary functions for which the diffference of the two compositions exactly vanishes:
\begin{align} g_1(x)&=\log(1+x),\\ g_2(x) &= e^x-1.\end{align}
as then $$g_1\circ g_2 = \log e^x = x = e^{\log(1+x)}-1 = g_2\circ g_1.$$
Of course, any function-inverse pair will work, like $2\sin(x)$ and $\arcsin(\frac12 x)$. The factor of $2$ shows that you don't need to be close to the identity map itself as guessed in Yves' answer, but this of course still requires the existence of an inverse.
Analysis of a general case
Aside from these perfect pairs,
the following basic calculation gives examples of arbitrary rates of vanishing at $x=0$.
Let $h_1,h_2$ be distinct (elementary, if you want) analytic functions around $0$ whose first $N$ taylor coefficients at $0$ vanish. Then let's consider $$g_1(x) = x+h_1(x), \\ g_2(x) = x+h_2(x).$$
For instance, $h_i$ could be $\tan$ and $\sin$ minus their respective $N$th Taylor polynomials. The tan-sin example you give is of this form with $N=2$ (since their $x^2$ coefficient is zero).
Their compositions are
$$ g_1(g_2(x)) = x+h_2(x) +h_1(x+h_2(x)),$$
$$ g_2(g_1(x)) = x+h_1(x) +h_2(x+h_1(x)).$$
and the subtraction is
$$\Delta:=g_1\circ g_2(x) - g_2\circ g_1(x) = \bigg(  h_1(x+h_2(x))- h_1(x) \bigg) - \bigg(h_2(x+h_1(x))- h_2(x)\bigg). $$
Taylor's theorem gives on a sufficiently small ball around $0$,
$$ h_i(x+h_j(x)) = h_i(x)+h_i'(x)h_j(x)+O(h_j(x)^2),$$
so
$$\Delta =h_1'(x)h_2(x) - h_2'(x)h_1(x) + O(h_1(x)^2+h_2(x)^2).$$
Since these terms vanish to high order （specifically, to order $O(x^{2N+1})$ the vanishing to high order is proven.
Optimal degree of vanishing
For the tan-sin example, this only proves $O(x^5)$, which is worse than the observed $O(x^7)$. In general, the order of vanishing can be improved by a small calculation to $O(x^6)$ (i.e. $O(x^{2N+2})$ in general) but no better. If you further write
$$ h_i(x) = a_i x^{N+1}+ O(x^{N+2})$$
then observe that $$ h_1'(x)h_2(x) = (N+1)a_1 a_2 x^{2N+1} + O(x^{2N+2}) $$
which has the same lowest order term as $h_2'(x)h_1(x)$. So these perfectly cancel, leaving you with an $O(x^{2N+2})$ error term as promised.
This is sharp in the class of functions I considered above, as can be seen by the example
$$ g_1(x)=e^x-1-\frac{x^2}2, \quad g_2(x) = \sin x + x^4, \quad \Delta = \frac{25}{144}x^6 + O(x^7).$$
Proving the $O(x^7)$ bound
Just a little more elbow grease explains the $O(x^7)$ term.
If you used the higher order Taylor expansion
$$ h_i(x+h_j(x)) = h_i(x)+h_i'(x)h_j(x)+h_i''(x)h_j(x)^2/2+O(h_j(x)^3),$$
we see that the terms $h_i''(x)h_j(x)^2/2+O(h_j(x)^3)$ only contribute to high order error terms $O(x^{2N+3})$. So the only term that could be only $O(x^{2N+2})$ is the term $h_i'(x)h_j(x)$. Expanding this to one more term than before:
$$ h_i(x) = a_i x^{N+1}+ b_ix^{N+2} +O(x^{N+3}),$$
we see
$$h_1'h_2(x)= (N+1)a_1a_2 x^{2N+1} + ((N+1)a_1b_2 + (N+2)a_2b_1)x^{2N+2} + O(x^{2N+3}).$$
Note that the second coefficient $(N+1)a_1b_2 + (N+2)a_2b_1$ is not cancelled when we subtract $h_2'(x)h_1(x)$, in contrast with the lower order term. So disregarding miracles like the perfect inverse pairs, we basically need $b_1,b_2=0$ to remove the  $x^{2N+2}$ term. This is precisely the case for the sin-tan example, as (recall $N=2$) there are no $bx^{N+2}=bx^4$ terms in the expansions of $\sin$ and $\tan$.
A: This occurs because the two functions are very similar to  $x$ for small $x$, so that the given composition cancels them efficiently.
Consider
$$\sin x\approx\arctan x\approx x,$$
then
$$\sin(\arctan x)-\arctan(\sin x)$$
is flat. The Taylor expansion is also of order $7$.
https://www.wolframalpha.com/input/?i=Taylor+expansion+of+sin%28arctan+x%29-arctan%28sin+x%29

For a nonic case, combine
$$\frac x{1-x^2}$$ and $$x-x^3+2x^5.$$
A: For any two functions $f(x)=x+a_1x^3+a_2x^5+a_3x^7+O(x^9)$, $g(x)=x+b_1x^3+b_2x^5+b_3x^7+O(x^9)$, we have
$$\begin{align}
f(g(x))=x&+(a_1+b_1)x^3+(3a_1b_1+a_2+b_2)x^3\\
&+(3a_1b_1^2+3a_1b_2+5a_2b_1+a_3+b_3)x^7+O(x^9),\\
g(f(x))=x&+(a_1+b_1)x^3+(3a_1b_1+a_2+b_2)x^3\\
&+(3a_1^2b_1+5a_1b_2+3a_2b_1+a_3+b_3)x^7+O(x^9).
\end{align}$$
