It's a generic omnidirectional motion model which isn't all that complicated, but it's a bit much to write out here. In this paper, it corresponds to the robot type (3,0), but the model in robot_localization is in 3D. Each variable has a nonlinear transition function, and we encode those transition functions as a matrix. So, for example, the transition function for x is x_new = x_old + vx * t + 0.5 * a * t^2, but v and a are given in the body frame of the robot, so we have to apply 3D rotations to them. The orientations have similar logic. Velocities (for x) are just vx_new = vx_old + at. Each row in the matrix is a transition function. We then derive the Jacobian analytically.
Naturally, this won't necessarily match ever robot perfectly, but in practice, it works pretty well for a broad range of robots. You an always increase the process_noise_covariance matrix values if the model doesn't match your robot well.