Need to understand Variational Gaussian Process in gpflow

Question

When I look at source of Base Likelihood, I saw this:

def variational_expectations(
        self, X: TensorType, Fmu: TensorType, Fvar: TensorType, Y: TensorType
    ) -> tf.Tensor:
        r"""
        Compute the expected log density of the data, given a Gaussian
        distribution for the function values,

        i.e. if
            q(f) = N(Fmu, Fvar)

        and this object represents

            p(y|f)

        then this method computes

           ∫ log(p(y=Y|f)) q(f) df.

        This only works if the broadcasting dimension of the statistics of q(f) (mean and variance)
        are broadcastable with that of the data Y.

        :param X: input tensor
        :param Fmu: mean function evaluation tensor
        :param Fvar: variance of function evaluation tensor
        :param Y: observation tensor
        :returns: expected log density of the data given q(F)
        """
        return self._variational_expectations(X, Fmu, Fvar, Y)

AFAIK, q(f)=N(f|μ,Σ) and Σ is covariance matrix, but gpflow define q(f) = N(Fmu, Fvar) when Fvar's shape is (N,1). What's that mean? Can you explain this to me? Thanks a lot!