Rotate objects around their origin along a global axis (scripted) without bpy.ops
TL;DR > Skip to last paragraph

Each object has it's own World Matrix. It's a 4x4 transform matrix that stores the object's final location, rotation and scale. By doing math operations directly on this matrix we can transform the object how ever we want.

Every world matrix can be decomposed into it's components. We will get a


location vector (size 3)
rotation quaternion (size 4)
scale vector (size 3)

loc, rot, scale = obj.matrix_world.decompose()



Now about rotations. A rotation is either represented by a (3x3 or 4x4) Rotation Matrix (Euler or Matrix class in Blender), Quaternion vector (Quaternion class) or an axis (Vector of size 3) with rotation value (radians or degrees). They can be all converted between each other, but we just need the 4x4 rotation matrix. There are different options to get it:

from mathutils import Matrix, Euler, Quaternion

# directly
Matrix.Rotation(angle, 4, axis)
# converting from Euler matrix
Euler((angleX, angleY, angleZ), 'XYZ').to_4x4()
# converting from quaternion
Quaternion((w,x,y,z)).to_matrix().to_4x4()


Transformations are done with matrix multiplication and they can be stacked. The order in which they are multiplied is important. So:

matA * matB * matC


means first we will transform with matC, then matB and last matA (it's reversed!). All the transformations are applied in global space. Here you can see what it looks like to do Translation * Rotation and Rotation * Translation:



The first end-result we can substitute with a local/global translation and then local rotation, and the second we can substitute with a local/global rotation and then local translation.

This means if the matrices are rotation matrices, only the first rotation matrix (and last applied) will result in global space rotation (so matA), because all the others are influenced by it.

This is how a World Matrix is composed. The order is again important:

matrix_world = matLoc * matRot * matScale


So to alter it and add an extra global rotation, we need to sneak a rotation matrix before matRot (so it's applied last and in global space) and after matLoc:

import bpy
from math import radians
from mathutils import Matrix

# we will demonstrate on an active object
obj = bpy.context.active_object

# define the rotation
rot_mat = Matrix.Rotation(radians(angle_in_degrees), 4, 'X')   # you can also use as axis Y,Z or a custom vector like (x,y,z)

# decompose world_matrix's components, and from them assemble 4x4 matrices
orig_loc, orig_rot, orig_scale = obj.matrix_world.decompose()
orig_loc_mat = Matrix.Translation(orig_loc)
orig_rot_mat = orig_rot.to_matrix().to_4x4()
orig_scale_mat = Matrix.Scale(orig_scale[0],4,(1,0,0)) * Matrix.Scale(orig_scale[1],4,(0,1,0)) * Matrix.Scale(orig_scale[2],4,(0,0,1))

# assemble the new matrix
obj.matrix_world = orig_loc_mat * rot_mat * orig_rot_mat * orig_scale_mat