Bezier curvature extrema For a planar cubic Bezier curve $B (x(t),y(t))$, I would like to find the values of parameter $t$ where the curvature (or curvature radius) is greatest/smallest. The formula for curvature is:
$$r = \dfrac{(x'^2+y'^2)^{(3/2)}}{x' (t) y''(t) - y'(t) x''(t)}$$
The problem is that there is that square root in it so I was wondering whether it is not possible to express the curvature extrema by some combination of the curve derivatives?
The idea is that for finding, say, the values where the slope of the curve is parallel to x-axis one needs to solve the quadratic function of the curve's first derivative. So I thought that maybe the extremes are in fact the values where the acceleration along the curve (or something similar) is greatest/smallest and I hoped that it would be possible to write it down as a polynomial.
 A: Short summary
The points of extremal curvature which are not inflection points will correspond to roots of the quintic polynomial
$$(x'''y' - x'y''')(x'^2 + y'^2)-3(x'x'' + y'y'')(x''y' - x'y'')$$
Finding these roots generally requires numeric methods; there can be no simple formula using radicals (i.e. fractional exponents) nor a ruler and compass construction of the points in question.
Avoiding square roots

The problem is that there is that square root in it

For this I'd indeed go with John's comment: the extrema of the curvature are extrema of the squared curvature. The converse may not hold if the curvature passes through zero, but if you are interested in the extrema of the absolute value of the curvature, then all maxima of the squared curvature will be relevant. So look at
$$k^2=\frac{(x'\,y'' - y'\,x'')^2}{(x'^2+y'^2)^3}$$
Characterizing extrema
Compute the derivative of that with respect to $t$, concentrate on the numerator of the resulting rational function, and you will get a polynomial of degree $7$ in $t$. It contains two factors: one of degree $2$ which is equivalent to $x'\,y'' - y'\,x''$ and describes the inflection points where the curvature is zero. And one of degree $5$ which gives you the remaining points of extremal curvature.
A concrete and illustrated example
I tried to determine whether all of these points are indeed relevant, and found several situations where all seven critical values for $t$ were real numbers between $0$ and $1$. So far I haven't found a really beautiful setup, where all seven points are easy to see, with sufficient distance between them and sufficient differences in curvature. Here is the best I have so far:
$$P_1=(0,4),P_2=(5,4),P_3=(7,0),P_4=(5,3)\\
t\in\{0.211, 0.326, 0.508, 0.634, 0.734, 0.852, 0.920\}$$
The curves looks like depicted below. Control points in cyan, inflection points in magenta and other extremal curvature points in red. I also included the osculating circles (light red) resp. tangent lines (light magenta) at the critical points, but they are all so close together that it's really hard to see.

Below is what the curvature looks like for the example above. It clearly depicts the fact that the curvature hardly differs at all between some of these extremal points, while at one maximum it's so large that it's way off the graph. I'm sorry I don't have a more balanced example.

Simple closed form impossible
Since there are at least five (and depending on your definition perhaps even seven) relevant solutions, I doubt you can get around solving this numerically. And if you do tackle this numerically, the above approach is probably easiest to understand and implement, so I'd not try looking for a different description. But more on that in the next section.
To be more precise: the five non-inflection points of extremal curvature in the example above are the roots of
$$18056\,t^5 - 48720\,t^4 + 49540\,t^3 - 23550\,t^2 + 5200\,t - 425 = 0$$
The Galois group of this polynomial is $S_5$, so there can be no way to describe its roots using radicals. Since most ruler and compass geometric constructions can't perform operations more complex than taking square roots, that also rules out reasonably simple geometric descriptions of these points.
Numeric approach
According to Wikipedia (and thanks to the comment by Jean Marie), the derivatives of your point can be computed like this:
\begin{align*}
B(t)&=(1-t)^3\,P_1+3(1-t)^2t\,P_2+3(1-t)t^2\,P_3+t^3\,P_4\\
B'(t)&=3(1-t)^2\,(P_2-P_1)+6(1-t)t\,(P_3-P_2)+3t^2\,(P_4-P_3)\\
B''(t)&=6(1-t)\,(P_3-2P_2+P_1)+6t\,(P_4-2P_3+P_2)\\
B'''(t)&=6(P_4-3P_3+3P_2-P_1)\\
B''''(t)&=0
\end{align*}
According to my computer algebra system, and unless I made a mistake, the derivative of the squared curvature is
$$\frac{\mathrm d\,k^2}{\mathrm dt}=
2\frac{\bigl(x''y' - x'y''\bigr)\bigl((x'''y' - x'y''')(x'^2 + y'^2)
-3(x'x'' + y'y'')(x''y' - x'y'')\bigr)}{(x'^2 + y'^2)^4}$$
Since you want this to be zero, you can ignore the denominator here. If you don't care about the inflection points, you can omit the first parenthesis in the numerator as well. If you know the coordinates of your four points, you can compute the coefficients of $t$ in the individual polynomials, combine them to the resulting polynomial (of degree 5 if you dropped the inflection condition or 7 if not) and then find roots of that polynomial numerically.
It is interesting to note that $(x'''y' - x'y''')$ has degree $1$ in $t$ and $(x''y' - x'y'')$ has degree $2$ in $t$, even though the degrees of the involved terms would suggest a higher resulting degree. But the terms in the highest coefficient cancel out, otherwise the second big parenthesis in the numerator would be a polynomial of degree $6$. If you do this in integer arithmetic, you will be able to cancel a common factor of $162$ from all the coefficients of this polynomial, but for floating point arithmetic it doesn't really matter.
Responses to comments
I'll try to answer some questions from your comment.

why at most 7 roots?

Because the relevant condition can be expressed as a polynomial of degree $7$ so there can never be more than these.

If I image a cubic Bezier I have a hard time finding seven places where the curvature could be most extreme.

I also have a hard time finding one where all the different points and curvatures can be seen well. But the example above should illustrate the kind of situation this is about. Whether you consider five or seven solutions depends on whether you care about inflection points with zero curvature or not.

Also, perhaps there are other ways to get the positions than squaring the curvature formula?

There may well be. But the solutions have to be equivalent. At least in the example above, we know that there are five resp. seven relevant points, so other solution cannot end up with a simpler description of these solutions. At least not simpler in terms of the degree of the describing polynomial.

And lastly, would something like the dot product of first and second derivative work? If not, why? That would result polynomial of degree 3.

Well, that's simply due to the definition of the curvature. While it sounds plausible that this dot product of yours would share some properties of curvature, it's a different quantity with different extremal parameters, so it does not fit the established concept of curvature as the inverse radius of osculating circles.
If you can't use that dot product all by itself, you could instead use it as part of a larger formula. Towards that goeal, the determinant (or 2d cross product, or whatever you want to call it) is proportional to the sine of the angle, just as the dot product is proportional to the cosine. So generally speaking you have
$$\langle v,w\rangle^2+(v\times w)^2=\lVert v\rVert^2\cdot\lVert w\rVert^2$$
for arbitrary vectors $v$ and $w$. Which means that in the case at hand you could write the squared curvature as
$$k^2=\frac{(x'^2+y'^2)(x''^2+y''^2)-(x'x''+y'y'')^2}{(x'^2+y'^2)^3}$$
but you gain nothing by doing this. The actual increase in degree comes from computing the derivative of such an expression, in order to find the extrema.
