I had a quick look at the bag file you posted (thanks for doing that). One thing I noticed is that you are not fusing any sources of absolute yaw orientation. In your first EKF (odom one), you do actually specify that the yaw should be fused from your imu, but by setting the imu_differential parameter to true, robot_localization will only look at the difference between yaw measurements. So what happens is the r_l nodes initialize the filter estimate with an orientation of unit quaternion (i.e., 0 deg yaw), and all subsequent imu data which impacts yaw is only changing it relative to the starting estimate of 0 degrees.