Are Riemannian submersions affine? Let $f:M \to N$ be a Riemannian submersion between two Riemannian manifolds $M,N$.
Consider $df $ as a section of the vector bundle $T^*M \otimes f^*TN$. Endow this bundle with the tensor product connection obtained from the Levi-civita connections on $TM$, $TN$ (formally, we first pullback the connection from $TN$ to $f^*TN$, then tensoring). 
Is $df$ a parallel section of $T^*M \otimes f^*TN$? ($\nabla df=0$?)
If not, is it true that $df^Tdf$ is a parallel section of $T^*M \otimes TM$?
(Of course, $df$ is parallel implies $df^Tdf$ is parallel but not vice versa)
 A: No, $df$ is not necessarily a parallel section of $T^{*}M \otimes f^{*}(TN)$ and $df^T \circ df$ is not a parallel section of $T^{*}M \otimes TM$. To see a concrete example, let us consider the projection from a surface of revolution onto a meridian. Explicitly, let $\gamma \colon I \rightarrow \mathbb{R}^3$ be a curve of the form $\gamma(v) = (\varphi(v), 0, \psi(v))$ with $\varphi(v) > 0$ and $\| \dot{\gamma}(v) \| = 1$. Let $X \colon I \times \mathbb{R} \rightarrow \mathbb{R}^3$ be the map
$$ X(v,\theta) := (\varphi(v) \cos(\theta), \varphi(v) \sin(\theta), \psi(v)). $$
The image of $X$ is the surface of revolution $S$ obtained by revolving $\gamma$ around the $z$-axis and $X$ is an immersion onto $S$. Pulling back the standard metric from $\mathbb{R}^3$, the induced metric on $M := I \times \mathbb{R}$ is $g = dv^2 + \varphi(v)^2 d\theta^2$. The projection onto the first coordinate $f(v,\theta) = v$ is then a Riemannian submersion when the metric on the interval $N := I$ is the standard metric. 
The non-vanishing covariant derivatives using Levi-Civita connection with respect to the frame $\partial_{v}, \partial_{\theta}$ are readily computed to be
$$ \nabla_{\partial_v} \partial_{\theta} = \nabla_{\partial_{\theta}} \partial_v = \frac{\varphi'(v)}{\varphi(v)} \partial_{\theta},\\ \nabla_{\partial_{\theta}} \partial_{\theta} = -\varphi(v) \varphi'(v) \partial_{v}. $$
Let us denote the coordinate on the interval $I$ by $x$ and the Levi-Civita connection on $I$ by $\nabla'$. Then $df = dv \otimes f^{*}(\partial_x)$ and $\nabla' \partial_x = 0$. Finally, denote by $\tilde{\nabla}$ the induced connection from $\nabla,\nabla',f$ on $T^{*}M \otimes f^{*}(TN)$. Computing, we have
$$ \tilde{\nabla}_{\partial_v}(df)(\partial_v) = f^{*}(\nabla')_{\partial_v}(df(\partial_v)) - df(\nabla_{\partial_v} \partial_v) = \nabla'_{\partial_x} \partial_x \circ f = 0, \\
\tilde{\nabla}_{\partial_v}(df)(\partial_{\theta}) = f^{*}(\nabla')_{\partial_v}(df(\partial_{\theta})) - df(\nabla_{\partial_v} \partial_{\theta}) = -df \left( \nabla_{\partial_v} \partial_{\theta} \right) =0,\\
\tilde{\nabla}_{\partial_{\theta}}(df)(\partial_{v}) = f^{*}(\nabla')_{\partial_{\theta}}(df(\partial_{v})) - df(\nabla_{\partial_{\theta}} \partial_{v}) = \nabla'_{df(\partial_{\theta})} \partial_x \circ f = 0,\\
\tilde{\nabla}_{\partial_{\theta}}(df)(\partial_{\theta}) = f^{*}(\nabla')_{\partial_{\theta}}(df(\partial_{\theta})) - df(\nabla_{\partial_{\theta}} \partial_{\theta}) = \varphi(v) \varphi'(v) f^{*}(\partial_x) $$
so
$$ \tilde{\nabla}(df) = \varphi(v) \varphi'(v) d\theta \otimes d\theta \otimes f^{*}(\partial_x). $$
The covariant derivative $\tilde{\nabla}(df)$ vanishes at $(v_0,\theta)$ if and only if $\varphi'(v_0) = 0$ if and only if the parallel $v = v_0$ is a geodesic and then it vanishes along this parallel. This occurs for example if $\varphi'(v) \equiv 0$ (in which case our surface is the flat cylinder) but doesn't occur generically.
Similarly, one can verify that $df^T \circ df = dv \otimes \partial_v$ and compute
$$ \tilde{\nabla}_{\partial_{\theta}} (df^T \circ df) = \frac{\varphi'(v)
}{\varphi(v)} dv \otimes \partial_{\theta} + \varphi(v) \varphi'(v) d\theta \otimes \partial_{v} $$
and this also doesn't vanish. Note that the expression above is a self-adjoint operator (with respect to $g$) because it is represented with respect to the orthonormal frame $\left( \partial_v, \frac{\partial_{\theta}}{\varphi(v)} \right)$ by the symmetric matrix
$$ \begin{pmatrix} 0 & \varphi'(v) \\ \varphi'(v) & 0 \end{pmatrix} $$
as expected from the fact that this is a derivative of the positive operator $df^T \circ df$.

BTW, a very similar computation can be done for a general Riemannian submersion. Given a Riemannian submersion $f$, O'Neill defined two tensors $T$ and $A$ which measure in some sense how non-trivial the submersion is. The tensor $T$ vanishes identically if and only if the fibers of $f$ are totally geodesic submanifolds of $M$. The tensor $A$ vanishes identically if and only if the horizontal distribution is integrable. In general, one can express $\tilde{\nabla} df$ using $df$ and the tensors $A$ and $T$. If they both vanish, then $df$ will be a flat section but not in general. In the calculation of $\tilde{\nabla} df$ above, the two middle terms actually vanish because the tensor $A$ vanishes identically (since the horizontal distribution is one-dimensional and so automatically integrable) and the first term (which calculates $\tilde{\nabla} df$ on two horizontal vector fields) will always vanish. The fourth term (which calculates $\tilde{\nabla} df$ on two vertical vector fields) will vanish if $T = 0$.
A: I don't think so.
Consider $M = N = \mathbb R^2$ equipped with the flat metric. Let $\theta : \mathbb R^2 \to \mathbb R$ be some smooth function we'll precise later, and consider the map
$$
df :=
\begin{pmatrix}
\cos \theta & \sin \theta
\\
-\sin \theta & \cos \theta
\end{pmatrix}.
$$
This is an orthogonal map for all $x,y$ by construction. I claim that we can find a smooth map $f : \mathbb R^2 \to \mathbb R^2$ such that its differential is this $df$; we can just pick a primitive of each factor separately, for example $f_{11}(x,y) = \int_0^x \cos \theta \, dx$, as there are no constraints between the factors.
For the right kind of $\theta$, the image of the map $f$ will contain an open set (we're trying to avoid stuff like $\theta(x,y) = x+y$ where the image is a circle). Between the preimage of this open set and itself, $f$ will be a map with discrete fibers, such that $df$ is an isometry. But since all the metrics involved are flat, $df$ being parallel means that $df$ should be constant, which we constructed it not to be.
