Using Parks-Mcclellan (or others) for FIR differentiator? I have never seen the use of Parks-Mcclellan for differentiator - most uses tend to be lowpass filters. Is this because Parks-Mcclellan implementation of FIR differentiator has numerical problems and cannot be successful?

It seems though that Parks-Mcclellan for a differentiator should be more studied and obvious, because in terms of fit of a function, a differentiator has frequency response that is basically a sloped straight line or an identity function.

In general, I am having some difficulty finding established optimal ways to implement FIR differentiators. Are there other known ways that are commonly used for FIR differentiators? What about lowpass differentiators, where frequency response basically is a differentiator up to some frequency and from some frequency to $\pi$ it is filtered out, with transition band in the middle?
There is no fundamental difference between the design of standard frequency selective filters (such as low pass, high pass, etc.) and differentiators. The only difference is that for a (linear phase) differentiator, the desired frequency response must be formulated as a purely imaginary function of frequency, unlike standard frequency selective filters, for which the desired response is purely real-valued. But that's a minor issue, and all important algorithmic details remain the same. Generally, algorithms only become different as soon as the desired response must be defined by a complex-valued function (with non-zero real and imaginary parts).

E.g., the Matlab function firpm, which implements the Parks McClellan method, can design FIR differentiators if you supply differentiator as filter type. The same is true for the least squares filter design function firls.

In sum, designing linear phase FIR differentiators is no more difficult than designing any other FIR filter.

Concerning low pass differentiators, take a look at this answer.