Relation between original points and 1st/2nd derivative points
Numerical signals are discrete and finite. So you are likely to hit issues either on end-points or discrete locations. I won't treat end-points here.

If you use only one operator, many signal processing folks care relatively little about a shift, as long as it is known, and constant. Here, the two-point derivative yields a $0.5$-point shift. You can consider that the derivative is valid at some mid-point. And you can interpolate it an integer locations if needed. Smoothing first can help limit a too-high fluctuations of differences.

A classical option is to choose an odd-length derivative operator. The three-point derivative writes:

$$ d_3[k] = (p[k+1]-p[k-1])/2$$

and gives you an information at the center of the schemes. Indeed, it coincides with the average of the left- and the right two-points skewed differences:  $ d_2[k] = p[k]-p[k-1]$ and $ d_2[k+1] = p[k+1]-p[k]$, since:

$$d_3[k] = (d_2[k] +d_2[k+1])/2$$

Many domains uses such schemes, and name them differently. To better find resources, some uses other terms, like generic ones, like "Finite difference coefficients", or specific ones, like "five-point stencil" or $n$-point gradient/Laplacian. One example is:

$$ d_5[k] = (p[k+2]- 8p[k+1] + 8p[k-1]-p[k-2])/12$$

The above thinking extends to higher order derivatives. One can further constraints methods using data properties, knowledge of noise, additional penalties, etc. To start with,  some literature on numerical derivatives:


CSERD Numerical Differentiation
Finite difference coefficients