I foolishly assumed that the OP was using 14.0 and didn't even run the original code. After plotting the function and seeing that it was real only for -1 < c < +1, I (wrongly) assumed that did the trick.
So to get to the point: My answer is not necessary and yours works for both 14.0.0.0 and 13.3.1.0.
For Mathematica 12.3.1.0 adding in a restriction works:
Subscript[M, r] = Maximize[{F[r, c], -0.9539 < c < 0.9539}, c]
(* {0.993421, {c -> 0.698707}} *)
Also for 14.0.0.0 you can get an "exact" answer (in terms of a Root object) if you rationalize the constants,
$a$ and
$r$.
Jim
