rho_buf = np.concatenate((tmp1, np.concatenate(
            (tmp2, rho, tmp2), axis=0), tmp1), axis=1)