How to create a 2D array in Python where each column has a different and random number of elements

I'm writing a Monte Carlo simulation of photons arriving at a detector. The total number of detected photons N_detected follows a Poisson distribution (which I'm generating using scipy.stats.poisson) and the time of detection of every single photon follows a given probability distribution function (PDF). In order to generate the time of detection of the N_detected photons I generate N_detected random numbers between 0 and 1 with the numpy.random.random() function and use the Inverse Transform Method.

I have to run the simulation several times. So, in order to avoid iterating a lot of times, I'd like to do every simulation at once using numpy arrays. As a final result, I'd like to obtain a 2D array of N_sim columns where every column corresponds to a different simulation and contains the generated times of detection. The problem is that each simulation produces a different number of photons (since it's random) and I don't find the way to create the 2D arrray with columns of different lengths.

One idea I had is creating all the columns with the same length (the largest N_deteceted) and fill the not needed elements with NaN and the rest with the random numbers I need, but I don't know I could do that.

This is the best I've been able to do so far:

import numpy as np
import numpy.ma as ma
from scipy.stats import poisson

beta = 0.0048    # Parameter of the PDF

""" Inverse of the time CDF (needed for the Inverse transform method) """

def inverse_time_cdf(x):
    return -np.log( (np.exp(-1000*beta)-1)*x + 1 )/beta


""" Simulation of N_sim experiments through the inverse transfrom method """

T = 1000    # Duration of the experiment [s]
rate = 3.945    # [events/s]
lam = rate*T    # Poisson distribution parameter

def photon_arrival_simulation(N_sim):
    N_detection = poisson.rvs(lam, size = N_sim)    # Number of detections in each experiment
    t = np.array([inverse_time_cdf(np.random.random(N)) for N in N_detection])

    return N_detection, t

If possible, I wanted to avoid the loop used in the list comprehension of the photon_arrival_simulation() function and also obtain a 2D array instead of an array of arrays (since I can't do array operation such as taking a log on an array of arrays).

I don't know if I should be posting this question here or in Physics Stack Exchange, but thanks in advance to anyone.