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 qubit-model QNodes, as well as all Gaussian continuous-variable 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.
Finite-difference: In the few cases where analytic gradients are not supported, for instance, if the QNode contains a non-Gaussian gate (e.g., a Kerr interaction or cubic phase gate) in-between the gate to be differentiated and the expectation value {}^1, then PennyLane falls-back to the finite-difference 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 non-Gaussian gate precedes the gate that is differentiated, then analytic differentiation is still supported.