Efficient multidimensional numerical integration in R and C++

Question

I'm trying to perform a 4-dimensional numerical integration in R using a function I wrote in C++ code which is then sourced in R using the Rcpp package.

Below there is my code:

// [[Rcpp::depends(RcppEigen)]]
// [[Rcpp::depends(RcppNumerical)]]
#define EIGEN_PERMANENTLY_DISABLE_STUPID_WARNINGS
#include <Eigen/Eigen>
#include <Rcpp.h>
#include <math.h>
#include <iostream>
#include <cstdlib>
#include <boost/math/special_functions/bessel.hpp>
#include <boost/math/special_functions/gamma.hpp>

using namespace std;
using namespace Rcpp;
using namespace Numer;

// [[Rcpp::depends(BH)]]

// [[Rcpp::export]]
double gaussian_free_diffusion(double x,double x0, double sigma, double t) {
  double pi = 2 * acos(0.0);
  double a1 = (1/sqrt(2.0 * pi * sigma * t));
  double b1 = exp(-pow((x - x0), 2.0)/(2.0 * sigma * t));
  double res = a1 * b1;
  return res;
}

// [[Rcpp::export]]
double integral_CppFunctionCUBATURE(NumericVector Zt_pos, const double& xA0, const double& xB0, const double& yA0,
   const double& yB0, const double& t1, const double& sigma){

  double xAt_pos = Zt_pos[0];
  double xBt_pos = Zt_pos[1];
  double yAt_pos = Zt_pos[2];
  double yBt_pos = Zt_pos[3];

  double temp_pbxA = gaussian_free_diffusion(xAt_pos, xA0, sigma,t1);
  double temp_pbxB = gaussian_free_diffusion(xBt_pos, xB0, sigma, t1);
  double temp_pbyA = gaussian_free_diffusion(yAt_pos, yA0, sigma,t1);
  double temp_pbyB = gaussian_free_diffusion(yBt_pos, yB0, sigma, t1);

  return (temp_pbxB * temp_pbyB) * (temp_pbxA * temp_pbyA);

};

integral_CppFunctionCUBATURE is the function I use with cubintegrate. In R, I would then do:

sourceCpp('./myCppcode.cpp')

xA0<-3
xB0<-2
yA0<-5
yB0<-10
t<-500
sigma<-1


cubature::cubintegrate(integral_CppFunctionCUBATURE,lower=rep(-1000,4), upper=rep(1000,4),
                       xA0=xA0, xB0=xB0, yA0=yA0,yB0=yB0, t1=t, sigma=sigma,method='cuhre')$integral

Result: 0.9999978

The multidimensional integration takes about 4 sec, but I would like to speed it up as much as possible, since when I run the multidimensional integration on several xA0 values this can take quite some time. Do you have any suggestion on how I could improve the speed of the code?

And also, would you suggest an alternative way to perform fast 4-dimensional integration in C++/Rcpp?