I’m also throwing some ideas in the pot, but you may be interested in using the particle trajectories functionality (although this isn’t exactly covering what you’re interested into). Otherwise, I’ve had some chance with a pseudo-code as follows (to be adapted to your use case)

ds1 = yt.load("/path/to/output_78")
ds2 = yt.load("/path/to/other/output")

old_ids = ds1.r["PartType0", "particle_identifier"]
old_order = np.argsort(old_ids)

def _old_density(field, data):
    ids = data["PartType0", "particle_identifier"]

    # find ids in old_ids
    inds = old_order[np.searchsorted(old_ids, ids, sorter=old_order)]

    # this is just to make sure I haven't messed up the ids matching
    np.testing.assert_allclose(ids, old_ids[inds])

    # Return the value of the density for those matched indices
    return ds1.r["density"][inds]

ds2.add_field(("PartType0", "old_density"), function=_old_density, sampling_type="cell", units=…)