Generate random directed connected graph using networkx?

Question

Context

While trying to generate a directed random connected graph in networkx, I am experiencing some difficulties combining the randomness with the "connected" property.

I use the following parameters to specify the properties of the graph:

1. size [2,100]: the number of nodes of the graph.
2. density (0,1]: The probability that an arbitrary edge is created.

Issue/Doubts

I am in doubt about the term connected in the context of a directed graph. What I mean is that any arbitrary node in the graph is able to reach all other nodes by traveling in the direction of the edges. I believe this does not mean that every edge between 2 nodes has to be bi-directional, as a detour may be possible to still reach nodes on the other side of an opposing-one-way-edge.

Another issue is that below certain edge densities, directed connected graphs are not possible anymore. However, I am not yet sure how I can account for that in the random graph generation. Additionally, for some densities, a directed graph may be possible, yet I would not directly know which edges to "randomly" select such that it actually is connected. Hence, I was wondering whether there exists an algorithm, or built in networkx function for this purpose.

Question

How can I generate a random graph of size=size that is directional and connected with an edge density=density in networkx?

MWE

The following MWE generates random graphs of a certain size and density, however I think they do not satisfy the connected property:

def gnp_random_connected_graph(
    density,
    size,
    test_scope,
):
    """Generates a random undirected graph, similarly to an Erdős-Rényi graph,
    but enforcing that the resulting graph is conneted.

    :param density:
    :param size:
    :param test_scope:
    """
    random.seed(test_scope.seed)
    edges = combinations(range(size), 2)
    G = nx.DiGraph()
    G.add_nodes_from(range(size))
    if density <= 0:
        return G
    if density >= 1:
        return nx.complete_graph(size, create_using=G)
    for _, node_edges in groupby(edges, key=lambda x: x[0]):
        node_edges = list(node_edges)

        random_edge = random.choice(node_edges)  # nosec - using a random seed.
        G.add_edge(*random_edge)
        for e in node_edges:
            if random.random() < density:  # nosec - no security application.
                G.add_edge(*e)