Monotonic, symmetrical soft-clipping polynomial
Okay, so to satisfy requirements 5, 6, 7, and 8, the soft-clipping function will have the form:

$$f(x) =
\begin{cases}
-1 &amp; x \le -1 \\
K \int\limits_{0}^{x} \big(1 - u^2 \big)^N \ du \quad &amp; -1\le x\le 1 \\
+1 &amp; 1 \le x \\
\end{cases}$$

where $K$ is a constant judiciously chosen to satisfy requirement 4.

The 1st derivative of $f(x)$ is

$$ f'(x) = K \big(1 - x^2 \big)^N $$

The 2nd derivative of $f(x)$ is

$$ f''(x) = K N \big(1 - x^2 \big)^{N-1} (-2x)$$

The 3rd derivative of $f(x)$ is

$$\begin{align}
 f'''(x) &amp;= K N(N-1) \big(1 - x^2 \big)^{N-2} (-2x)^2 \ + \ K N \big(1 - x^2 \big)^{N-1} (-2) \\
&amp;= K \big(1 - x^2 \big)^{N-2} \bigg( N(N-1)(-2x)^2 - 2N \big(1 - x^2 \big)  \bigg) \\
\end{align}$$

and, for $n \ge 1$, the $n$th derivative is

$$ f^{(n)}(x) = K \big(1 - x^2 \big)^{N-n+1} \, g_n(x) $$

where $g_n(x)$ is some ($n$-1)th order polynomial function of $x$ and is finite in value.  This can be proven inductively by considering the ($n$+1)th derivative:

$$\begin{align}
\\
f^{(n+1)}(x) &amp;= \tfrac{d}{dx} \Big( f^{(n)}(x) \Big) \\
&amp;= \tfrac{d}{dx} \Big(K \big(1 - x^2 \big)^{N-n+1} \, g_n(x)\Big) \\
&amp;= K(N-n+1)\big(1 - x^2 \big)^{N-n} (-2x) g_n(x) \, + \, K \big(1-x^2 \big)^{N-n+1}g'_n(x) \\
&amp;= K\big(1 - x^2 \big)^{N-n} \Big( (N-n+1)(-2x) g_n(x) \, + \, (1-x^2) g'_n(x) \Big) \\
&amp;= K\big(1 - x^2 \big)^{N-n} \, g_{n+1}(x) \\
\\
\end{align}$$

where $ g_{n+1}(x) = (N-n+1)(-2x) g_n(x) + (1-x^2) g'_n(x)  $.

Because of differentiation, the order of polynomial $g'_n(x)$ is one less than the order of $g_n(x)$, but polynomial $(1-x^2)g'_n(x)$ is one order greater than $g_n(x)$ and so also is $(-2x) g_n(x)$.

When $x = \pm 1$, then the first $N$ derivatives are zero, 
$$ f^{(n)}(x) \Bigg|_{x=\pm 1} = 0 \qquad \text{for } 1 \le n \le N $$
making this polynomial maximally flat at $x = \pm 1$.

The integrand is a binomial and can be expressed as a power series using binomial expansion:

$$\begin{align}
\big(1 - u^2 \big)^N &amp; = \sum\limits_{n=0}^{N} \binom{N}{n} \big(-u^2\big)^n (1)^{N-n} \\
  &amp; = \sum\limits_{n=0}^{N} \frac{N!}{n!(N-n)!} \big(-u^2\big)^n (1)^{N-n} \\
  &amp; = \sum\limits_{n=0}^{N} \frac{N!}{n!(N-n)!} (-1)^n u^{2n} \\
\end{align}$$

So the integral can be expressed as an integral of a power series:

$$\begin{align}
\int\limits_{0}^{x} \big(1 - u^2 \big)^N \ du &amp; = \int\limits_{0}^{x} \sum\limits_{n=0}^{N} \frac{N!}{n!(N-n)!} (-1)^n u^{2n} \ du \\
&amp; = \sum\limits_{n=0}^{N} \frac{N!}{n!(N-n)!} (-1)^n \int\limits_{0}^{x} u^{2n} \ du \\
&amp; = \sum\limits_{n=0}^{N} \frac{N!}{n!(N-n)!} \frac{(-1)^n}{2n+1} u^{2n+1} \Bigg|_0^x \\
&amp; = \sum\limits_{n=0}^{N} \frac{N!}{n!(N-n)!} \frac{(-1)^n}{2n+1} x^{2n+1} \\
&amp; = x \sum\limits_{n=0}^{N} \frac{N!}{n!(N-n)!} \frac{(-1)^n}{2n+1} \big(x^2 \big)^n \\
\end{align}$$

When $x = \pm 1$, we get

$$\begin{align}
 \int\limits_{0}^{\pm 1} \big(1 - u^2 \big)^N \ du &amp;= \pm \sum\limits_{n=0}^{N} \frac{N!}{n!(N-n)!} \frac{(-1)^n}{2n+1} \\
  &amp;= \pm \frac{4^N}{2N+1}\binom{2N}{N}^{-1}\\
\end{align}$$

(there is a little improperness in the integral limit.)  The bottom equality is because of this Appendix proof.  So the scaler $K$ must be

$$\begin{align}
 K &amp;= \left( \frac{4^N}{2N+1}\binom{2N}{N}^{-1} \right)^{-1} \\
\\
 &amp;= \frac{2N+1}{4^N}\binom{2N}{N} \\
\\
 &amp;= \frac{2N+1}{4^N} \frac{(2N)!}{(N!)^2} \\
\\
 &amp;= \frac{(2N+1)!}{4^N (N!)^2} \\
\end{align}$$

This makes the entire soft-clipping function to be:

$$ f(x) = \begin{cases} 
-1 &amp; x \le -1 \\
\\
x \cdot \frac{(2N+1)!}{4^N \ N!} \sum\limits_{n=0}^{N} \frac{(-1)^n}{(2n+1)\,n!\,(N-n)!} \big(x^2 \big)^n \quad &amp; -1 \le x \le +1 \\
\\
+1 &amp; +1 \le x \\
\end{cases} $$

and it appears that the odd-power (that is the power = $2n+1$) polynomial coefficients are

$$ a_n = (-1)^n \frac{(2N+1)!}{4^N \, N! \,(2n+1) \, n! \, (N-n)!} $$

The polynomials (without splicing to the $\operatorname{sgn}(x) = \pm 1$ saturated components) look like



I think the order $2N+1$ starts at 1 and goes to 9 (or $0 \le N \le 4$)

With the saturation attached, the curves look like



The soft-clipping function is continuous everywhere and all derivatives, up to the $N$th derivative are continuous everywhere and the ($N$+1)th derivative and higher are continuous everywhere except at the splices where $x = \pm 1$.

Here are the same set of curves but with the scaling adjusted so that the slope around $x=0$ (or "gain") remains at 1 (or "0 dB gain").  



The point of discontinuity (where the polynomial is spliced to a constant $\pm 1$) is at $x= \pm K$.  I really don't think one needs to get over a $7$th-order ($N=3$) softclipper.  Oversampling by a factor $4$x suffices to deal with this $7$th-order polynomial and prevent any aliasing to the original baseband before upsampling.