The negative connected correlation function because of the momentum conservation

cr_ideal.py

import numpy as np
from scipy.spatial.distance import pdist
import matplotlib.pyplot as plt
from scipy.stats import binned_statistic
from numba import njit


@njit
def pairwise_dot(x):
    """
    Calculate the pair-wise dot product of quantity x. This function should
        be called for calculating the spatial correlation of quantity x.
        (~3x faster with numba on my machine)

    Args:
        x (:obj:`numpy.ndarray`): the quantity for pairwise dot calculation
            the shape of x is (N, dimension)

    Return:
        :obj:`numpy.ndarray`: the pair-wise dot product of x as a 1D array.
            The order of the result is the same as the pairwise distance
            calculated from :obj:`scipy.spatial.distance.pdist`, or :obj:`pdist_pbc`
    """
    N = len(x)
    result = np.empty(int((N * N - N) / 2))
    idx = 0
    for i in range(N):
        for j in range(i+1, N):
            result[idx] = x[i] @ x[j]
            idx += 1
    return result


def random_sample(N, L=10, spd=1, frames=1000, D=3, want_flctn=True):
    """
    Calculate the correlation function of the velocities of ideal gas.

    The positions and velocities of the gas were randomly sampled in each frame.

    All particles have the same speed, instead of following the MB distribution.

    flctn = fluctuation

    Args:
        N (:obj:`int`): the number of particles in the system.
        L (:obj:`float`): the length of the size of the PBC.
        spd (:obj:`float`): the speed of all particles in different frames.
        frames (:obj:`int`): the total simulation frame.
        D (:obj:`int`): the dimension of the system.
        want_flctn (:obj:`bool`): if True calculate the connected correlation
            otherwise calculate the disconnected correlation

    Return:
        :obj:`tuple`: all the pairwise distances and
            the pairwise dot product of velocity fluctuations
    """
    distances = []
    velocity_flctn = []
    for f in range(frames):
        # draw the positions
        positions = np.random.uniform(0, L, (N, D))
        # draw the velocities
        velocities = np.random.normal(0, 1, (N, D))
        velocities /= np.linalg.norm(velocities, axis=1)[:, np.newaxis]
        # calculate the fluctuation of the velocities
        if want_flctn:
            flctn = velocities - velocities.mean(axis=0)[np.newaxis, :]
        else:
            flctn = velocities
        distances.append(pdist(positions))
        velocity_flctn.append(pairwise_dot(flctn))
    return np.concatenate(distances), np.concatenate(velocity_flctn)


if __name__ == "__main__":
    np.random.seed(0)
    N = 50
    frames = 1000
    # calculate & plot the connected correlation frunction of velocity
    dists, velocities = random_sample(N, frames=frames, want_flctn=True)
    bins = np.linspace(0, 5, 26)
    bc = bins[:-1]
    mean, _, _ = binned_statistic(dists, velocities, bins=bins, statistic="mean")
    std, _, _ = binned_statistic(dists, velocities, bins=bins, statistic="std")
    count, _, _ = binned_statistic(dists, velocities, bins=bins, statistic="count")
    plt.errorbar(bc, mean, std / np.sqrt(count), marker='o', capsize=2, label='connected', color='tomato')
    # calculate & plot the disconnected correlation frunction of velocity
    dists, velocities = random_sample(N, frames=frames, want_flctn=False)
    bins = np.linspace(0, 5, 26)
    bc = bins[:-1]
    mean, _, _ = binned_statistic(dists, velocities, bins=bins, statistic="mean")
    std, _, _ = binned_statistic(dists, velocities, bins=bins, statistic="std")
    count, _, _ = binned_statistic(dists, velocities, bins=bins, statistic="count")
    plt.errorbar(bc, mean, std / np.sqrt(count), marker='o', capsize=2, label='disconnected', color='teal')
    plt.plot([0, bc[-1]], [0, 0], color='k')
    plt.xlabel('r', fontsize=16)
    plt.ylabel('C(r)', fontsize=16)
    plt.legend(fontsize=16)
    plt.tight_layout()
    plt.show()

Here is the result for N=50