$\newcommand{\al}{\alpha}
\newcommand{\de}{\delta}
\newcommand{\De}{\Delta}
\newcommand{\ep}{\varepsilon}
\newcommand{\ga}{\gamma}
\newcommand{\Ga}{\Gamma}
\newcommand{\la}{\lambda}
\newcommand{\Si}{\Sigma}
\newcommand{\thh}{\theta}
\newcommand{\om}{\omega}
\newcommand{\R}{\mathbb{R}}
\newcommand{\Z}{\mathbb{Z}}
\newcommand{\F}{\mathcal{F}}
\newcommand{\E}{\operatorname{\mathsf E}}
\renewcommand{\P}{\operatorname{\mathsf P}}
\newcommand{\ii}[1]{\operatorname{\mathsf I}\{#1\}}
\newcommand{\tf}{\widetilde{f}}$

I will borrow the excellent idea of user fedja about a simplex. Without loss of generality (wlog), the iid random variables (rv's) are simple functions, of the form $\sum_{j=1}^k x_j\,1_{E_j}$ for some disjoint measurable sets $E_j$ and some (not necessarily distinct) $x_j$'s in $[0,1]$. Then $\E S$ is a convex function of $x:=(x_1,\dots,x_n)\in[0,1]^n$ and hence attains its maximum in $x\in[0,1]^n$ at some $x\in\{0,1\}^n$. Hence, wlog each of the iid rv's takes only two values, $1$ and $0$, with probabilities $p$ and $q:=1-p$, for some $p\in[0,1]$. So, the exact upper bound on $\E S$ is
\begin{equation*}
\E S_*:=\sup_{p\in[0,1]}\E|N_1(p)-N_2(p)|=\sup_{p\in[0,1]}T(p),\tag{0}
\end{equation*}
where $N_1(p),N_2(p)$ are iid binomial rv's with parameters $n,p$,
\begin{equation*}
T(p):=
\sum_{i,j=0}^n|i-j|\binom ni\binom nj a_{i+j}(p),
\end{equation*}
and
\begin{equation*}
a_k(p):=a_{n,k}(p):=p^{k}q^{2n-k}.
\end{equation*}

We have
\begin{equation*}
T(p)=\sum_{k=0}^{2n}a_k(p)s_k,
\end{equation*}
\begin{align*}
s_k:=\sum_{i=0}^k|i-(k-i)|\binom ni\binom n{k-i}
&= \frac2n\,m \binom{n}{m} (n-k+m) \binom{n}{k-m} \tag{1} \\
&= \left\{
\begin{aligned}
2n\,\binom{n-1}{m-1}^2 &\text{ if }k=2m-1,\\
2n\,\binom{n-1}{m-1} \binom{n-1}{m} &\text{ if }k=2m,
\end{aligned}
\right.
\end{align*}
where $m:=\left\lceil \frac{k}{2}\right\rceil$; the second equality in (1) was obtained with Mathematica (I will check it later).

Next,
\begin{equation*}
T(p)=T_1(p)+T_0(p), \tag{1.5}
\end{equation*}
\begin{equation*}
T_1(p):=\sum_{m=1}^{n}a_{2m-1}(p)s_{2m-1},\quad
T_0(p):=\sum_{m=1}^{n}a_{2m}(p)s_{2m}.
\end{equation*}

Further, $T_1(p)$ is a hypergeometric expression and, accordingly, it satisfies the diff. eq.
\begin{equation*}
pq (1-4npq) T_1 '(p)+(p-q) (1-2npq) T_1 (p)
+(p-q)(pq)^2 T_1 ''(p)=0. \tag{2}
\end{equation*}
Also, we have the symmetry $T_1(1-p)=T_1(p)$. So, $T'_1(1/2)=0$.
Moreover, if $0<p<1/2$, $2npq>1$, and $T_1 '(p)$, then (2) and the positivity of $T_1(p)$ imply $T_1 ''(p)>0$. So, there are no local maxima of $T_1$ in the interval
\begin{equation*}
J_n:=\{p\colon 0<p<1/2,2npq>1\}.
\end{equation*}

Let us now state two lemmas.

Lemma 1. For $n\ge12$ and $p\in[0,1/2)\setminus J_n$, we have $T_1(p)<T_1(1/2)$, so that
\begin{equation*}
\max_{0<p<1}T_1(p)=T_1(1/2). \tag{3}
\end{equation*}

Lemma 2. For $n\le11$, (3) holds.

These lemmas will be proved at the end of the answer.

Lemmas 1 and 2, together with the previous observations that $T_1$ is positive and symmetric and that there are no local maxima of $T_1$ on $J_n$, yield (3) for all natural $n$.

The consideration of $T_0(p)$ is similar to, and much easier than, that of $T_1(p)$. We have
\begin{equation*}
(1-4npq) T_0'(p)-2n(p-q) T_0(p)+ (p-q)pq T_0''(p)=0.
\end{equation*}
So, there are no local maxima of $T_0$ in the entire interval $[0,1/2)$. Also, $T_0(0)=0<T_0(1/2)$. Therefore, $\max_{0<p<1}T_0(p)=T_0(1/2)$, and so, in view of (3), $\max_{0<p<1}T(p)=T(1/2)$, for all natural $n$. That is, the exact upper bound on $\E S$ is
\begin{equation*}
\E S_*=T(1/2)=\E|N_1(1/2)-N_2(1/2)|.
\end{equation*}

For large $n$, using the normal approximation to the binomial distribution with parameters $n,p=1/2$, we see that the best upper bound on $\E S$ can be approximated as follows:
\begin{equation*}
\E|N_1(1/2)-N_2(1/2)|
\approx\E|Z_1-Z_2|\sqrt{n\tfrac12\,\tfrac12}=\E|Z\sqrt2|\sqrt{n/2}=\frac1{\sqrt\pi}\,\sqrt n,
\end{equation*}
where $Z,Z_1,Z_2$ are iid standard normal rv's.

It remains to prove the two lemmas. Lemma 2 is proved by direct calculation, using the fact that $T_1(p)$ is a polynomial in $p$ of degree $2n$. It remains to present

*Proof of Lemma 1.* We have
\begin{equation*}
T_1(p)\le T(p)=\E|N_1(p)-N_2(p)|\le\sqrt{\E(N_1(p)-N_2(p))^2}=
\sqrt{2\mathsf{Var}\,N_1(p)}=\sqrt{2npq}\le1 \tag{4}
\end{equation*}
for $p\in[0,1/2)\setminus J_n$. On the other hand,
\begin{equation*}
T_1(1/2)=
2^{1-2 n} n\sum_{m=1}^{n}\binom{n-1}{m-1}^2
=2^{1-2 n} n \binom{2 (n-1)}{n-1},
\end{equation*}
which is easily seen to be increasing in natural $n$, with $T_1(1/2)=1.0091\ldots>1$ for $n=12$.
So, in view of (4), Lemma 1 follows.