def photon_number_mean(mu, cov, j, hbar=2):
    """ Calculate the mean photon number of mode j of a Gaussian state.

    Args:
        mu (array): vector of means of the Gaussian state using the ordering
            :math:`[q_1, q_2, \\dots, q_n, p_1, p_2, \\dots, p_n]`
        cov (array): the covariance matrix of the Gaussian state
        j (int): the j :sup:`th` mode
        hbar (float): the ``hbar`` convention used in the commutation
            relation :math:`[q, p]=i\\hbar`

    Returns:
        float: the mean photon number in mode :math:`j`.
    """
    num_modes = len(mu) // 2
    return (mu[j] ** 2 + mu[j + num_modes] ** 2 + cov[j, j] + cov[j +
        num_modes, j + num_modes] - hbar) / (2 * hbar)
