Hi @kareem_essafty, that is correct, the gradients are computed by sampling the quantum hardware directly, as in the paper you linked
In general, PennyLane has two approaches to calculating the gradient on quantum hardware:

Analytic gradients: these are supported by all qubitmodel QNodes, as well as all Gaussian continuousvariable circuits. This is done by sampling the quantum devices at two additional parameter values, proving the exact gradient.
The derivation for the equations used to determine the analytic gradient are given in the paper linked above.

Finitedifference: In the few cases where analytic gradients are not supported, for instance, if the QNode contains a nonGaussian gate (e.g., a Kerr interaction or cubic phase gate) inbetween the gate to be differentiated and the expectation value {}^1, then PennyLane fallsback to the finitedifference approximation.
This is still computed via sampling the quantum device, but is a numerical approximation determined by sampling the circuit at infinitesimal shifts around the parameter.
{}^1 if the nonGaussian gate precedes the gate that is differentiated, then analytic differentiation is still supported.