Simulation of Wiener Field (Python)

Question

I am interested in simulating Wiener field in dimension 2 over [0,1]^{d}, see definition here :https://encyclopediaofmath.org/wiki/Wiener_field

I am working with Python.

I have done the naïvest thing i could think of and this somehow work but when i take a small step size for discretization this become really long to compute. The construction of the covariance matrix is awful in what i have done (the four imbricate "for" loops are kind of prohibitive). How can i improve it ?

n=100 # level of discretization direction-wise
# construct the covariance matrix
C=np.zeros((n,n,n,n))
for i in range(n):
  for j in range(n):
    for k in range(n):
      for l in range(n):
        C[i][j][k][l]=min(i,k)*min(j,l)/(n*n)

C2=np.reshape(C,(n*n,n*n))
#Generate the 2D Wiener field
x = np.random.multivariate_normal(np.zeros(n*n), C2, 1)
W=np.reshape(x, (n,n))