Because you are working in one dimension, there's a polynomial-time algorithm, using dynamic programming. Assume $x_1,\dots,x_n$ are sorted in increasing order. Notice that each cluster must contain a contiguous sequence of $x$'s, say $x_i,x_{i+1},\dots,x_j$. Let $f(n_0,k_0)$ denote the minimum of the maximum inter cluster distance for dividing the points $x_1,\dots,x_{n_0}$ into $k_0$ clusters. Then you can express $f(n_0,k_0)$ in terms of $f(m,k_0-1)$ for $m<n_0$:
$$f(n_0,k_0) = \min_m \max(f(m,k_0-1), \text{variance}(x_{m+1},x_{m+2},\dots,x_{n_0})),$$
where $m$ ranges over $m=1,2,\dots,n_0-1$.
Consequently, you can evaluate $f(\cdot,\cdot)$ at all $O(nk)$ inputs in $O(n^2k)$ time. You can improve this to $O(nk)$ time by setting $g(n_0,k_0)=\min(f(1,k_0),\dots,f(n_0,k_0))$ and noting
$$\begin{align*}
f(n_0,k_0) &= \max(g(m,k_0-1), \text{variance}(x_{m+1},x_{m+2},\dots,x_{n_0}))\\
g(n_0,k_0) &= \min(g(n_0,k_0-1), f(n_0,k_0)).
\end{align*}$$
You will need to use the fact that you can compute the mean and variance of $x_{m+1},\dots,x_{n_0}$ from the mean and variance of $x_1,\dots,x_m$ and the mean and variance of $x_1,\dots,x_{n_0}$. You can compute the mean and variance of $x_1,\dots,x_{n_0}$ for all $n_0$ in $O(n)$ time by maintaining a running sum and running sum of squares.