Is there a way to read the rotation_quaternion data of an object that is under the control of the physics engine?
There's a fairly simple way to get the transformations of an object under a physics simulation without baking (which also applies to other situations such as parent-child relationships, constraints, etc).

To do this you simply need to access the object's world matrix property. Here are some examples you can test through the interactive python console:

# Accessing through the location property will give the initial value before the simulation starts
&gt;&gt;&gt; C.object.location
Vector((0.0, 0.0, 14.47101879119873))

# The world matrix though gives the real, simulated location
&gt;&gt;&gt; C.object.matrix_world.to_translation()
Vector((0.5996798276901245, 1.5040934085845947, 1.7373782396316528))

# Same goes for rotation
&gt;&gt;&gt; [ degrees(a) for a in C.object.rotation_euler ]
[0.0, 0.0, 0.0]

&gt;&gt;&gt; [ degrees(a) for a in C.object.matrix_world.to_euler() ]
[176.1907675277052, -11.192468775907107, -74.38630617504491]

# To access quaternion rotation values, simply use "to_quaternion"
&gt;&gt;&gt; C.object.matrix_world.to_quaternion()
Quaternion((0.17901940643787384, -0.7067779302597046, 0.6757702231407166, -0.10840550810098648))


You can also add an app handler that calculates and prints the simulated location and rotation after every frame change:

import bpy
from math import degrees

def print_simulated_loc_rot( scene ):
    C = bpy.context
    o = C.object

    simulated_loc   = o.matrix_world.to_translation()
    simulated_rot   = o.matrix_world.to_euler()
    sim_rot_degrees = [ degrees(v) for v in simulated_rot ]

    print_friendly_loc = [ round( v, 2 ) for v in simulated_loc   ]
    print_friendly_rot = [ round( v, 2 ) for v in sim_rot_degrees ]

    print( "Location: ", print_friendly_loc )
    print( "Rotation: ", print_friendly_rot )

bpy.app.handlers.frame_change_post.append( print_simulated_loc_rot )