The reason you haven't read about "rotation about an axis" in 4D is that there is no such thing. Perhaps I should say "the reason is that any rotation that leaves fixed a single axis also leaves fixed another axis orthogonal to it, and hence leaves fixed an entire plane." You were probably hoping to find something with the property that
$$
Rv = v
$$
for any $v$ in the "axis" $L$, and the property that
$$
Rw \cdot w = \cos \theta
$$
for any vector $w$ orthogonal to $L$, i.e., all vectors $w$ in the plane orthogonal to $L$ are rotated by an angle of $\theta$. But, alas, there is no such thing.
The trick is to understand that our standard way of describing rotations in 3-space is slightly flawed: instead of saying that "we rotated 20 degrees about the $z$ axis", for instance, we could say "we rotated in the $xy$-plane by $20$ degrees" (i.e., naming the plane of rotation rather than the single axis that happens to be orthogonal to the plane). If we did the latter, then generalizing would work fine: we could talk (in 4-space) about rotating in the $xy$-plane, which leaves two perpendicular dimensions ($z$ and $w$) fixed. Such rotations are relatively easy to write down, but they do not describe all possible rotations in 4-space: there are also ones that rotate by some amount in a plane spanned by vectors $v_1, v_2$, and by some (possibly different) amount in a plane spanned by vectors $u_1, u_2$, with the $u$s orthogonal to the $v$s.
It's not too tough to write down proofs of these claims, at least if you know about eigenvalues and eigenvectors; if you want, I can expand my answer to do so. But I think the main contribution that I can make is to help you understand that you're looking for something that (provably) does not exist, and perhaps you should attempt to gain a different understanding.
Post-comment additions:
@Widawensen asked that I expand a bit on my remarks, so I'll do so, giving proofs for a few things I've said.
First: Throughout, $R$ will be a rotation in $4$-space, or the $4 \times 4$ matrix of this rotation with respect to the standard basis.
Claim 1: If, for some nonzero vector $v$, we have $Rv = v$, then there's another vector $w$ orthogonal to $v$ with $Rw = w$. In other words: if a rotation fixes all points on a single axis (the span of the vector $v$), then it actually fixes all points on an entire plane.
Proof: A rotation matrix preserves lengths, so its eigenvalues must all be complex (or real numbers) of absolute-value 1. We know, since $Rv = v$ and $v$ is nonzero, that one of the eigenvalues is 1. So the eigenvalues are $1, a, b, c$. Either all of $a,b, c$ are real (in which case they're all $\pm 1$), or some of them are complex. In the real case, since the determinant of a rotation is $+1$, we must have $abc = +1$, hence either 0 or two are $-1$.
Case 1: all four eigenvalues are $+1$: then $R$ is the identity, and any vector orthogonal to $v$ can be chosen as $w$ and we're done.
Case 2: The evalues are $+1, +1, -1, -1$. Since $+1$ is an eigenvalue twice, the eigenspace $E$ corresponding to $+1$ is two-dimensional$^{**}$, and it evidently contains $v$. Extending $v$ to an orthogonal basis $
\{v, w\}$ of $E$ via Gram-Schmidt, we get our vector $w$ and we're done.
If there are complex eigenvalues, the situation is similar to the second case: Because the characteristic polynomial is a real polynomial (no imaginary numbers!) the complex eigenvalues occur in complex conjugate pairs, so we have
- Case 3: The evalues are $1, a, p + iq, p - iq$, where $p^2 + q^2 = 1$ because all eigenvalues must have length 1 for a rotation. The determinant is then $1 \cdot a \cdot (p^2 + q^2) = a$. But we know the det of a rotation is $1$, so $a = 1$. Hence the eigenspace for $1$ is again 2-dimensional, and the analysis for Case 2 applies.
Claim 2: Rotations that fix all points of some plane in 4-space do not actually exhaust all rotations of 4-space: there are rotations of 4-space that move the points of every 2D plane.
Proof: Look at the matrix
\begin{bmatrix}
-1 & 0 & 0 & 0 \\
0 & -1 & 0 & 0 \\
0 & 0 & c & -s \\
0 & 0 & s & c
\end{bmatrix}
where $c = s = \cos \frac{\pi}{4} = \frac{\sqrt{2}}{2}$.
The eigenvalues of this matrix are $-1, -1, c+is, c-is$. In particular, the number $1$ is not an eigenvalue, so there is no vector $v$ with $Rv = v$.
Claim 3: Every rotation $R$ in 4-space is a product $R = ST$, where $S$ is a rotation that fixes all points of a plane $A$, and $T$ is a rotation that fixes all points of a plane $B$, and the planes $A$ and $B$ are orthogonal, in the sense that for any vectors $a \in A$ and $b \in B$, we have $A \cdot b = 0$.
Proof: Again this is by cases, based on the eigenvalue structure. The five possible eigenvalue sets are $(1,1,1,1)$, $(1,1,-1, -1)$, $(-1, -1, -1, -1)$, $(1, 1, a+bi, a-bi)$ (where $b \ne 0$ and $a^2 + b^2 = 1)$, $(a+bi, a-bi, c+di, c-di)$, where $b, d \ne 0$ and $a^2 + b^2 = c^2 + d^2 = 1$. We'll examine these one at a time:
Case 1: $(+1, +1, +1, +1)$. $R$ is the identity, and can be factored as $I$ \cdot $I$, where the first $I$ fixes all points of the $xy$-plane, and the second fixes all points of the $zw$-plane, so that these two planes serve the roles of $A$ and $B$.
Case 2: $(+1, +1, -1, -1)$. The eigenspace for $+1$ is 2-dimensional, so pick an orthonormal basis $v_1, v_2$ for it. Extend this, by Gram Schmidt, to an orthonormal basis for 4-space, $v_1, v_2, v_3, v_4$.
Because the eigenspaces for $+1$ and $-1$ must be orthogonal, the span of $v_3, v_4$ must be the eigenspace for $-1$. Let $M$ be the matrix whose columns are $v_1, v_2, v_3, v_4$, and $D$ be the diagonal matrix with diagonal entries $-1, -1, +1, +1$. Then $S = M D M^t$ negates $v_1$ and $v_2$ and fixes $v_3$ and $v_4$. Picking $T = I$, we have $R = S T$, and we're done.
Case 3: $(+1, +1, a+bi, a-bi)$. The eigenspace for $+1$ is 2-dimensional, so pick an orthonormal basis $v_1, v_2$ for it. Extend this, by Gram Schmidt, to an orthonormal basis for 4-space, $v_1, v_2, v_3, v_4$. Let $M$ be the matrix whose columns are $v_1, v_2, v_3, v_4$. Look at
$$
Q = M R M^t
$$
This matrix fixes $e_1$ and $e_2$ and has the same eigenvalues as $R$ (because $M$ is a rotation). Since it's a rotation, its columns must be orthogonal. That means that $Q$ must have the form
$$
\begin{bmatrix}
1 & 0 & 0 & 0 \\
0 & 1 & 0 & 0 \\
0 & 0 & * & * \\
0 & 0 & * & *
\end{bmatrix}
$$
where the $2 \times 2$ block at the lower right is actually an 2D rotation. That means it has the form
\begin{bmatrix}
\cos \theta & -\sin \theta \\
\sin \theta & \cos \theta
\end{bmatrix}
for some $\theta$ and has eigenvalues $\cos \theta \pm i ~ \sin \theta$. But we know that its eigenvalues must be $a \pm bi$. So the actual form of $Q$ must be
$$
Q = \begin{bmatrix}
1 & 0 & 0 & 0 \\
0 & 1 & 0 & 0 \\
0 & 0 & a & -b \\
0 & 0 & b & a
\end{bmatrix}
$$
(or its transpose).
This matrix factors into $I \cdot Q$, where $I$ fixes the $zw$-plane, and where $Q$ fixes all points of the $xy$ plane. So we have
$$
M R M^t = IQ
$$
Conjugating by $M$, we get
\begin{align}
M^t (M R M^t) M &= M^t (IQ) M \\
(M^t M) R (M^t M) &= M^t (IMM^tQ) M \\
R &= (M^t IM) (M^tQ M)
\end{align}
which shows that $R$ is a product of something that fixes the $v_3, v_4$ plane (namely $M^t I M = I$) and something that fixes the $v_1, v_2$ plane (namely $M^t Q M = R$). The key point is that $M^t M = M M^t = I$, because its columns forma an orthonormal basis for 4-space.
Case 3: $(-1, -1, -1, -1)$. In this case, $R$ is just the negative identity, and factors nicely as
$$
R = \begin{bmatrix}
1 & 0 & 0 & 0 \\
0 & 1 & 0 & 0 \\
0 & 0 & -1 & 0 \\
0 & 0 & 0 & -1
\end{bmatrix}
\begin{bmatrix}
-1 & 0 & 0 & 0 \\
0 & -1 & 0 & 0 \\
0 & 0 & 1 & 0 \\
0 & 0 & 0 & 1
\end{bmatrix}.
$$
Case 4: Like case 2. $R$ is (up to conjugation by a matrix $M$ whose columns are its eigenvectors) just the matrix of "Claim 2" above.
Case 5: Similar, except that there's no single plane that's either pointwise fixed or pointwise negated: instead, there are two perpendicular planes that are each rotated by some non-multiple-of-$\pi$ amount.
Footnote: $^{**}$ FShrike points out that it's not obvious that just having $+1$ as an eigenvalue twice means that there are two eigenvectors with eigenvalue $1$. After all, consider the matrix $\pmatrix{1 & 1 \\ 0 & 1}$, right? In general, all we can say is that the generalized eigenspace for $\lambda = +1$ is two-dimensional. That gives us an eigenvector $v$ (which we might as well adjust to have length 1) for the eigenvalue $+1$, and a second independent vector $u$ for which (1) $Au \in span \{v, u\}$, (2) $(A-I)u$ may be nonzero, but (3) $(A-I)^2 u$ is zero. By applying Gram-Schmidt to the pair $\{v, u\}$, we get a pair $\{v, w\}$, where the three properties above also hold for $w$.
By item 1, we know that $Aw = cv + dw$ for some $c, d$, so $(A-I)w = cv + (d-1)w$. Applying this a second time gives us
\begin{align}
(A-I)^2 w
&= c (A-I)v + (d-a)(A-I)w \\
&= 0 + (d-1)[cv + (d-1)w]\\
&= c(d-1)v + (d-1)^2w
\end{align}
which must be $0$, and because $v$ and $w$ are independent, we see that $(d-1)^2 = 0$, i.e. $d = 1$. Now look at $Aw = cv + w$. Because $A$ is orthogonal (which means length-preserving), the left-hand side has length $\|Aw\| = \|w\| = 1$. The right-hand side (remember that $v$ and $w$ are orthogonal and unit length, because of Gram Schmidt) has length $\sqrt{c^2 + 1}$. Hence $c = 0$, and we've shown that the second generalized eigenvector, $w$, is actually an eigenvector, so we're done.
[As @FShrike notes, if you know the spectral theorem, this is all a good deal easier.]
Best Answer
Any symplectomorphism has Jacobian with determinant one. It is relatively easy to see that the Hamiltonian $H_n(x)=\Vert x\Vert^n$ for $n \geq 2$ has as flows the "rotations" you mention. (I.e., by powers of $r$.)
Indeed, $\nabla H_n(x)=n\Vert x\Vert^{n-2}x,$ thus letting $\Phi_t^n$ denote the Hamiltonian flow of $H_n$ at time $t$, we have that $\Phi^4_{1/4}$ is the rotation you mainly want in the question. For the generalized question of rotations by powers $r^n$ of $r$, just take $\Phi_{1/(n+2)}^{n+2}$.
This also generalizes for the case that Jyrki mentions, i.e. $\theta \circ r$. For this case, take $F$ to be an antiderivative for $f(r)=\theta(r) \cdot r$ and pick the Hamiltonian $H=F \circ \Vert \cdot \Vert$. Then $$\nabla H(x)=F'(\Vert x \Vert)\cdot \frac{x}{\Vert x \Vert}=\theta(\Vert x \Vert) \cdot \Vert x \Vert \cdot \frac{x}{\Vert x \Vert}=\theta(\Vert x \Vert) \cdot x.$$ Then the flow at time $1$ is the desired rotation.
It may not be clear why this is conceptual and natural without some background in symplectic geometry, so some words about that: