Improving performance of iterative 2D Numpy array with multivariate random generator

In a UxU periodic domain, I simulate the dynamics of a 2D array, with entries denoting x-y coordinates. At each time step, the "parent" entries are replaced by new coordinates selected from their normally distributed "offsprings", keeping the array size the same. To illustrate:

import numpy as np
import random
np.random.seed(13)

def main(time_step=10):
    def dispersal(self, litter_size_):
        return np.random.multivariate_normal([self[0], self[1]], [[sigma**2*1, 0], [0, 1*sigma**2]], litter_size_) % U

    U = 10
    sigma = 2.
    parent = np.random.random(size=(4,2))*U

    for t in range(time_step):
        offspring = []
        for parent_id in range(len(parent)):
             litter_size = np.random.randint(1,4)   # 1-3 offsprings reproduced per parent  
             offspring.append(dispersal2(parent[parent_id], litter_size))
        offspring = np.vstack(offspring)

        indices = np.arange(len(offspring))
        parent = offspring[np.random.choice(indices, 4, replace=False)]  # only 4 survives to parenthood

    return parent

However, the function can be inefficient to run, indicated by:

from timeit import timeit
timeit(main, number=10000)

that returns 40.13353896141052 secs.

A quick check with cProfile seems to identify Numpy's multivariate_normal function as a bottleneck.

Is there a more efficient way to implement this operation?