Simulate a 2D-beta distribution: spatial simulations

I'm doing spatial simulations and I would like to simulate a matrix that represent a truncated 2D-bimodal distribution with U shape. In other words, I would expect high density towards the borders of the matrix and lower density in the middle.

Looking around I have noticed that the distribution that better would fit is the beta distribution. However, do you have any idea of how I could reproduce a 2d beta u shaped density distribution?

I would like to simulate this matrix in R.