Gradient estimates Gradient estimates (and especially the differential Harnack) for harmonic functions on Riemannian manifolds were proved by Cheng and Yau in 1975, by using Bochner's formula. However, it seems that similar techniques had been used by Bernstein around 1910. My question is, what was the new contribution of Cheng-Yau from the PDE standpoint? Was it mainly the Riemannian aspect or was the technique new even for harmonic functions in a flat domain? In other words, did Bernstein derive the differential Harnack inequality by his method (It seems that he did derive the plain gradient estimates)?
 A: That part of Cheng-Yau's paper is a followup to Yau's paper "Harmonic functions on Riemannian manifolds" from the same year, adapting Yau's work on complete manifolds to geodesic balls. Yau's work is about the maximum principle applied to product-rule identities like
$$\Delta\frac{u+c_1}{\sqrt{|\nabla u|^2+c_2}}=\frac{\Delta u}{\sqrt{|\nabla u|^2+c_2}}+\cdots$$
So applicability of the maximum principle gives a pointwise estimate like $u+c_1\leq \kappa\sqrt{|\nabla u|^2+c_2}$ which is what is now called "differential Harnack" or (depending on context) Yau/Cheng-Yau/Li-Yau estimate. In Yau's paper, there is an application to a gradient estimate for harmonic functions on the unit ball in Euclidean space, by reinterpreting the solution as a solution of an elliptic equation on hyperbolic space. But I do not know if the estimate would be obvious by other means to an expert on elliptic PDE.
In the weaker integrated form, these estimates were well-known as Harnack inequalities, except that the Riemannian-geometric dependence of the constants (as pointed out in comment) was not understood. Moreover I'm not aware of any prior work where pointwise estimates like Yau's are achieved; at least, I have heard from experts that the directly analogous Li-Yau gradient inequality was new even for the standard heat equation $u_t=u_{xx}$. I believe Bernstein's gradient estimates are prototypes of the (uniform-type) Schauder estimates from 1930s and are not of this pointwise type, see e.g. "S. N. Bernshtein's contribution to the theory of partial differential equations" http://dx.doi.org/10.1070/RM1961v016n02ABEH004101. Such gradient estimates cannot be integrated to Harnack inequalities.
Yau/Cheng-Yau's technique of applying the laplacian to a composite expression has its origin in Bernstein's "method of auxiliary functions" to find a priori estimates. From this perspective some people might say that Yau/Cheng-Yau's contribution is only the discovery of a useful auxiliary function, but other people might call that view overly reductive. It would be very interesting to know what kinds of auxiliary functions were considered from (say) 1910-1960 but that is beyond my historical knowledge.
In Yau's paper, the setup for analysis is an extension of Omori's maximum principle to allow maximum principle analysis on complete spaces with  Ricci curvature bounded below. Cheng-Yau's paper gave another proof. The Omori-Yau principle has had several applications in geometry, but I don't know if it has consequences that can be advertised to analysts who don't care about geometry.
