All four functions call stats::acf() or stats::pacf() to do the calculations. They are simply wrappers to provide more convenient output or improved plots.

The help file for stats::acf() and stats::pacf() gives a brief explanation. For ACF:

For type = "correlation" and "covariance" , the estimates are based on the sample covariance. (The lag 0 autocorrelation is fixed at 1 by convention.)

and for PACF

The partial correlation coefficient is estimated by fitting autoregressive models of successively higher orders up to lag.max.

Later it gives the reference:

Venables, W. N. and Ripley, B. D. (2002) Modern Applied Statistics with S. Fourth Edition. Springer-Verlag.
(This contains the exact definitions used.)

If you look at the reference, you will find the standard covariance estimates of the ACF (p390), and the PACF defined using the Yule-Walker equations (on p398).

If you want to see precisely how it is done, you can always look at the code. After all, R is open source. The code for stats::acf() includes

acf <- .Call(C_acf, x, lag.max, type == "correlation")

which tells you that it calls C to do the calculation. If you look at the C sources, you will find the relevant function here:

Similarly, the PACF uses the C code here:

The last question looks like a homework problem. To solve moment problems for problems like this, backshift notation is helpful for re-arranging the equation into a suitable form -- you want y_t on the left and only lagged error terms on the right. Then you need a couple of standard results for infinite power series.

The reason I ask regarding autocorrelations is ESTIMA use Burg as default https://estima.com/webhelp/index.html?context=1290 see Correlate - method used to compute correlations - autocorrelations - Yule–Walker vs Burg Methods.

My second question, the process is: y(t) = mu + phi1*y(t-1) + beta*(t) + eps(t) , as from ADF test, third model form.