Langevin equation autocorrelation function

the average is over multiple Brownian particles, so experimentally, you would track these particles and compute the property you want at each time step for every particle and do the mean

The average here is over different realizations of random process $eta(tau)$. If we measured the values of noise only at discrete time instants $tau_1, tau_2, .., tau_n$, we could write a joint probability of these as $$ w(x_1, x_2, ..., x_n), $$ which in this case should be a multivariate Gaussian distribution with diagonal covariance matrix. As we take time interval to be smaller and smaller and pass to a continuous limit, we have to average over a functional $$ w[x(t)] $$ and use functional calculus. A good book covering this issues is the first volume by KLyatskin. (This is unfortunately a translation from an old Russian text. If somebody can recommend an equivalent text taht is mroe readily available, I would be glad toa dd it here.)

In the equation at hand the averaging is trivial, since the correlation function is known: $$ int_0^t dtauint_0^{t'} dtau'langle eta_i(tau)eta_j(tau')rangle e^{-gamma(t+t'-tau-tau')} = int_0^t dtauint_0^{t'} dtau'delta_{i,j}delta(tau-tau') e^{-gamma(t+t'-tau-tau')}= delta_{i,j}theta(t'-t)int_0^t dtau e^{-gamma(t+t'-2tau)} + delta_{i,j}theta(t-t')int_0^{t'} dtau' e^{-gamma(t+t'-2tau')}, $$ where $theta(t)$ is the Heaviside step function, which accounts for the overlap of integration ranges necessary for delta-function taking a non-zero value.