Creating a load carrying multirobot simulation in pybullet

Question

I am planning to do multirobot load carrying simulations which would like this where the two circles are robots and the rectangle is a load

I have found an excellent library that has a very simple holonomic point robot simulation https://github.com/maxspahn/gym_envs_urdf in which I have two robots loaded which have a urdf associated to them and are loaded into the environment.

To this simulation, I'm looking for ways in which I can add a load like the one shown in the figure.1, the load is attached to a fixed continuous joint which allows the robot to rotate in place with the load on the top.

To accomplish this would the urdf files have to modified or can this be done directly through the pybullet python API. Any help is deeply appreciated. Thank you.

code for the simulation, to run this, the library has to installed.

import gym
import numpy as np
from urdfenvs.robots.tiago import TiagoRobot
from urdfenvs.robots.generic_urdf import GenericUrdfReacher
from urdfenvs.robots.prius import Prius
# import sys
# sys.path.append("/home/josyula/Programs/MAS_Project/gym_envs_urdf/")
import pybullet
def run_multi_robot(n_steps=1000, render=False, obstacles=False, goal=False):
    robots = [
        GenericUrdfReacher(urdf="pointRobot.urdf", mode="vel"),
        GenericUrdfReacher(urdf="pointRobot.urdf", mode="vel"),
        # GenericUrdfReacher(urdf="ur5.urdf", mode="acc"),
        # GenericUrdfReacher(urdf="ur5.urdf", mode="acc"),
        # TiagoRobot(mode="vel"),
        # Prius(mode="vel")
    ]

    env = gym.make(
        "urdf-env-v0",
        dt=0.01, robots=robots, render=render
    )
    n = env.n()
    action = np.ones(n) * -0.2
    pos0 = np.zeros(n)
    pos0[1] = -0.0
    base_pos = np.array([
        [0.0, 1.0, 0.0],
        [0.0, -1.0, 0.0],
        [0.0, -2.0, 0.0]
        ])
    ob = env.reset(pos=pos0, base_pos=base_pos)
    print(f"Initial observation : {ob}")
    if goal:
        from examples.scene_objects.goal import dynamicGoal
        env.add_goal(dynamicGoal)

    # if obstacles:
    #     from examples.scene_objects.obstacles import dynamicSphereObst2
    #     env.add_obstacle(dynamicSphereObst2)
    pybullet.createVisualShape(pybullet.GEOM_BOX, halfExtents=[0.5, 0.5, 0.5])
    print("Starting episode")
    history = []
    for _ in range(n_steps):
        ob, _, _, _ = env.step(action)
        history.append(ob)
    env.close()
    return history

    #add load onto robots

if __name__ == "__main__":
    run_multi_robot(render=True, obstacles=True, goal=True)