diff --git a/crates/energy/src/conservation.rs b/crates/energy/src/conservation.rs
index 0b38620..6f3e6f3 100644
--- a/crates/energy/src/conservation.rs
+++ b/crates/energy/src/conservation.rs
@@ -1,4 +1,5 @@
 use bevy::prelude::*;
+use forces::prelude::{RotationalWorkEvent, WorkDoneEvent};
 
 /// Enum representing different types of energy
 #[derive(Debug, Clone, Copy, PartialEq, Eq, Hash, Component, Reflect)]
@@ -86,7 +87,7 @@ pub struct EnergyTransferEvent {
 /// Component for precise energy accounting
 #[derive(Component, Debug, Reflect)]
 #[reflect(Component)]
-pub struct EnergyAccountingLedger {
+pub struct EnergyBalance {
     /// History of all transactions, newest first
     pub transactions: Vec<EnergyTransaction>,
     /// Maximum number of transactions to store
@@ -116,7 +117,7 @@ pub struct EnergyTransaction {
     pub duration: f32,
 }
 
-impl Default for EnergyAccountingLedger {
+impl Default for EnergyBalance {
     fn default() -> Self {
         Self {
             transactions: Vec::new(),
@@ -127,7 +128,7 @@ impl Default for EnergyAccountingLedger {
     }
 }
 
-impl EnergyAccountingLedger {
+impl EnergyBalance {
     /// Record a new energy transaction
     pub fn record_transaction(&mut self, transaction: EnergyTransaction) {
         match transaction.transaction_type {
@@ -258,6 +259,72 @@ impl EnergyDriftMonitor {
     }
 }
 
+/// System to ensure entities with Mass have energy balance tracking
+pub fn initialize_energy_balance(
+    mut commands: Commands,
+    query: Query<Entity, (With<forces::prelude::Mass>, Without<EnergyBalance>)>,
+) {
+    for entity in query.iter() {
+        commands.entity(entity).insert(EnergyBalance::default());
+    }
+}
+
+/// System to track work done by forces and record in energy balance
+pub fn track_work_from_forces(
+    mut work_events: MessageReader<WorkDoneEvent>,
+    mut query: Query<&mut EnergyBalance>,
+    time: Res<Time>,
+) {
+    for event in work_events.read() {
+        if let Ok(mut ledger) = query.get_mut(event.entity) {
+            // Correct transaction direction based on work sign
+            let (transaction_type, source, destination) = if event.work > 0.0 {
+                (TransactionType::Input, None, Some(event.entity))
+            } else {
+                (TransactionType::Output, Some(event.entity), None)
+            };
+
+            ledger.record_transaction(EnergyTransaction {
+                transaction_type,
+                amount: event.work.abs(),
+                source,
+                destination,
+                timestamp: time.elapsed_secs(),
+                transfer_rate: event.work.abs() / time.delta_secs().max(f32::EPSILON),
+                duration: time.delta_secs(),
+            });
+        }
+    }
+}
+
+/// System to track rotational work done by torques and record in energy balance
+pub fn track_rotational_work_from_torques(
+    mut work_events: MessageReader<RotationalWorkEvent>,
+    mut query: Query<&mut EnergyBalance>,
+    time: Res<Time>,
+) {
+    for event in work_events.read() {
+        if let Ok(mut ledger) = query.get_mut(event.entity) {
+            // Same logic as linear work: positive work = energy input, negative = output
+            let (transaction_type, source, destination) = if event.work > 0.0 {
+                (TransactionType::Input, None, Some(event.entity))
+            } else {
+                (TransactionType::Output, Some(event.entity), None)
+            };
+
+            ledger.record_transaction(EnergyTransaction {
+                transaction_type,
+                amount: event.work.abs(),
+                source,
+                destination,
+                timestamp: time.elapsed_secs(),
+                transfer_rate: event.work.abs() / time.delta_secs().max(f32::EPSILON),
+                duration: time.delta_secs(),
+            });
+        }
+    }
+}
+
 /// Plugin to manage energy conservation systems
 pub struct EnergyConservationPlugin;
 
@@ -269,11 +336,22 @@ impl Plugin for EnergyConservationPlugin {
             .register_type::<EnergyQuantity>()
             .register_type::<TransactionType>()
             .register_type::<EnergyTransaction>()
-            .register_type::<EnergyAccountingLedger>()
+            .register_type::<EnergyBalance>()
             // Add resources
             .init_resource::<EnergyConservationTracker>()
             // Add event channel
-            .add_message::<EnergyTransferEvent>();
+            .add_message::<EnergyTransferEvent>()
+            // Add systems with explicit ordering
+            .add_systems(
+                Update,
+                (
+                    initialize_energy_balance,
+                    ApplyDeferred,
+                    (track_work_from_forces, track_rotational_work_from_torques)
+                        .after(forces::PhysicsSet::ApplyForces),
+                )
+                    .chain(),
+            );
     }
 }
 
@@ -284,7 +362,7 @@ mod tests {
     #[test]
     fn test_ledger_arithmetic() {
         // Test that net_energy_change = total_input - total_output
-        let mut ledger = EnergyAccountingLedger::default();
+        let mut ledger = EnergyBalance::default();
 
         ledger.record_transaction(EnergyTransaction {
             transaction_type: TransactionType::Input,
@@ -324,7 +402,7 @@ mod tests {
     #[test]
     fn test_current_flux_sums_active_rates() {
         // Test that current_flux sums transfer rates correctly
-        let mut ledger = EnergyAccountingLedger::default();
+        let mut ledger = EnergyBalance::default();
 
         let current_time = 10.0;
 
diff --git a/crates/energy/src/electromagnetism/charges.rs b/crates/energy/src/electromagnetism/charges.rs
index 34acd5f..6308762 100644
--- a/crates/energy/src/electromagnetism/charges.rs
+++ b/crates/energy/src/electromagnetism/charges.rs
@@ -197,10 +197,12 @@ pub fn apply_coulomb_pairwise_forces(
                 continue;
             }
 
-            // Coulomb force: F = k·q₁·q₂/r²
+            // Coulomb force: F_on_A = -k·q₁·q₂/r² · r̂_AB
+            // Same-sign charges (k_qq > 0) → repulsive (force pushes A away from B)
+            // Opposite-sign charges (k_qq < 0) → attractive (force pulls A toward B)
             let k_qq = config.coulomb_constant * charge_a * charge_b;
             let force_magnitude = k_qq / r.powi(2);
-            let force_bare = (force_magnitude / r) * r_vec; // F_bare = (k·q₁·q₂/r²) · r̂
+            let force_bare = -(force_magnitude / r) * r_vec;
 
             // Apply C¹ force-switch for smooth cutoff
             let switch = force_switch(r, config.switch_on_radius, config.cutoff_radius);
diff --git a/crates/energy/src/lib.rs b/crates/energy/src/lib.rs
index a661cb8..b4b0036 100644
--- a/crates/energy/src/lib.rs
+++ b/crates/energy/src/lib.rs
@@ -113,9 +113,9 @@ pub mod prelude {
     pub use super::{EnergySystem, EnergyTransferError};
 
     pub use crate::conservation::{
-        EnergyAccountingLedger, EnergyConservationPlugin, EnergyConservationTracker,
-        EnergyDriftMonitor, EnergyQuantity, EnergyTransaction, EnergyTransferEvent, EnergyType,
-        TransactionType, conversion_efficiency, verify_conservation,
+        EnergyBalance, EnergyConservationPlugin, EnergyConservationTracker, EnergyDriftMonitor,
+        EnergyQuantity, EnergyTransaction, EnergyTransferEvent, EnergyType, TransactionType,
+        conversion_efficiency, verify_conservation,
     };
 
     pub use crate::electromagnetism::prelude::*;
diff --git a/crates/energy/src/thermodynamics/thermal.rs b/crates/energy/src/thermodynamics/thermal.rs
index d7f0603..21b49a3 100644
--- a/crates/energy/src/thermodynamics/thermal.rs
+++ b/crates/energy/src/thermodynamics/thermal.rs
@@ -1,6 +1,4 @@
 use bevy::prelude::*;
-use forces::PhysicsSet;
-use forces::core::newton_laws::Mass;
 use matter::geometry::Radius;
 use std::collections::HashMap;
 use utils::{SpatialIndexSet, SpatiallyIndexed, UnifiedSpatialIndex, force_switch};
@@ -474,20 +472,22 @@ fn check_thermal_stability(
 /// **Conservation**: Thermal energy tracked, but not yet integrated with ledger.
 fn sync_thermal_energy(
     mut commands: Commands,
-    changed_temps: Query<(Entity, &Temperature, &HeatCapacity, &Mass), Changed<Temperature>>,
+    changed_temps: Query<(Entity, &Temperature, &HeatCapacity), Changed<Temperature>>,
 ) {
-    // **LP-0 SCAFFOLDING**: U = m·c_p·T (simplified thermal energy).
+    // **LP-0 SCAFFOLDING**: U = C·T (simplified thermal energy).
+    // HeatCapacity already includes mass: C = m·c_p (J/K).
     // Future: Use proper thermodynamic potentials (enthalpy for constant P, etc.).
 
-    for (entity, temp, heat_cap, mass) in changed_temps.iter() {
-        // Thermal energy: U = m·c_p·T (Joules)
+    for (entity, temp, heat_cap) in changed_temps.iter() {
+        // Thermal energy: U = C·T (Joules)
+        // where C = HeatCapacity = m·c_p (J/K), already accounts for mass.
         //
         // **APPROXIMATION**: Assumes constant c_p (valid for small ΔT).
         // Assumes T_ref = 0 K (absolute thermal energy).
         //
         // **IRL PHYSICS**: For phase changes, need enthalpy H = U + PV.
         // For proper thermodynamics, need U(S,V) or H(S,P).
-        let thermal_energy = mass.value * heat_cap.value * temp.value;
+        let thermal_energy = heat_cap.value * temp.value;
 
         commands.entity(entity).insert(EnergyQuantity {
             value: thermal_energy,
@@ -515,8 +515,18 @@ impl Plugin for ThermalSystemPlugin {
                 PreUpdate,
                 mark_thermal_entities_spatially_indexed.in_set(SpatialIndexSet::InjectMarkers),
             )
-            // Thermal transfer in Update - doesn't write AppliedForce, no conflict
-            .add_systems(Update, sync_thermal_energy);
+            // Thermal conduction → flush commands → sync energy.
+            // conduction inserts Temperature via Commands; apply_deferred flushes
+            // so sync_thermal_energy sees Changed<Temperature> in the same frame.
+            .add_systems(
+                Update,
+                (
+                    compute_fourier_conduction,
+                    ApplyDeferred,
+                    sync_thermal_energy,
+                )
+                    .chain(),
+            );
     }
 }
 
diff --git a/crates/energy/src/waves/propagation.rs b/crates/energy/src/waves/propagation.rs
index 50346ab..da480ac 100644
--- a/crates/energy/src/waves/propagation.rs
+++ b/crates/energy/src/waves/propagation.rs
@@ -6,9 +6,7 @@ pub(crate) struct WaveGrid(pub(crate) SpatialGrid);
 
 use super::normalize_or;
 use super::oscillation::{WaveParameters, angular_frequency, wave_number};
-use crate::conservation::{
-    EnergyAccountingLedger, EnergyQuantity, EnergyTransaction, TransactionType,
-};
+use crate::conservation::{EnergyBalance, EnergyQuantity, EnergyTransaction, TransactionType};
 
 // Calculate modified angular frequency with dispersion
 #[inline]
@@ -144,7 +142,7 @@ pub fn apply_wave_damping_with_energy(
         (Entity, &WaveParameters, Option<&mut EnergyQuantity>),
         With<WaveCenterMarker>,
     >,
-    mut ledger_query: Query<&mut EnergyAccountingLedger>,
+    mut ledger_query: Query<&mut EnergyBalance>,
 ) {
     let dt = time.delta_secs();
     if dt == 0.0 {
@@ -179,8 +177,8 @@ pub fn apply_wave_damping_with_energy(
                 }
             }
 
-            // Record transaction to ledger
-            if let Ok(mut ledger) = ledger_query.single_mut() {
+            // Record transaction to this entity's ledger
+            if let Ok(mut ledger) = ledger_query.get_mut(entity) {
                 ledger.record_transaction(EnergyTransaction {
                     transaction_type: TransactionType::Output,
                     amount: energy_lost,
diff --git a/crates/forces/src/core/gravity.rs b/crates/forces/src/core/gravity.rs
index 13ae3bd..1091fd4 100644
--- a/crates/forces/src/core/gravity.rs
+++ b/crates/forces/src/core/gravity.rs
@@ -437,15 +437,13 @@ pub fn calculate_barnes_hut_attraction(
 
     affected_query
         .par_iter_mut()
-        .for_each(|(entity, transform, _, mut force)| {
+        .for_each(|(entity, transform, mass, mut force)| {
             let position = transform.translation;
 
-            if bodies.iter().any(|&(e, _, _)| e == entity) {
-                return;
-            }
-
             let force_vector = calculate_barnes_hut_force(
+                entity,
                 position,
+                mass.value,
                 &quadtree.root,
                 theta,
                 softening,
@@ -457,7 +455,9 @@ pub fn calculate_barnes_hut_attraction(
 }
 
 pub fn calculate_barnes_hut_force(
+    affected_entity: Entity,
     affected_position: Vec3,
+    affected_mass: f32,
     node: &spatial::QuadtreeNode,
     theta: f32,
     softening: f32,
@@ -470,7 +470,8 @@ pub fn calculate_barnes_hut_force(
         let distance_squared = direction.length_squared();
         let softened_distance_squared = distance_squared + softening_squared;
         let force_magnitude =
-            gravitational_constant * node.mass_properties.total_mass / softened_distance_squared;
+            gravitational_constant * affected_mass * node.mass_properties.total_mass
+                / softened_distance_squared;
 
         if !force_magnitude.is_finite() {
             return Vec3::ZERO; // Skip non-finite BH aggregated force
@@ -482,7 +483,12 @@ pub fn calculate_barnes_hut_force(
     if node.children.iter().all(|c| c.is_none()) {
         let mut total_force = Vec3::ZERO;
 
-        for &(_, position, mass) in &node.bodies {
+        for &(entity, position, mass) in &node.bodies {
+            // Skip self-interaction
+            if entity == affected_entity {
+                continue;
+            }
+
             let direction = position - affected_position;
             let distance_squared = direction.length_squared();
 
@@ -491,7 +497,8 @@ pub fn calculate_barnes_hut_force(
             }
 
             let softened_distance_squared = distance_squared + softening_squared;
-            let force_magnitude = gravitational_constant * mass / softened_distance_squared;
+            let force_magnitude =
+                gravitational_constant * affected_mass * mass / softened_distance_squared;
 
             if !force_magnitude.is_finite() {
                 continue; // Skip non-finite BH leaf force
@@ -506,7 +513,9 @@ pub fn calculate_barnes_hut_force(
     let mut total_force = Vec3::ZERO;
     for child_node in node.children.iter().flatten() {
         total_force += calculate_barnes_hut_force(
+            affected_entity,
             affected_position,
+            affected_mass,
             child_node,
             theta,
             softening,
@@ -746,9 +755,16 @@ mod tests {
             params.barnes_hut_max_bodies_per_node,
         );
 
-        for (i, _) in entities.iter().enumerate() {
-            let bh_force =
-                calculate_barnes_hut_force(bodies[i].0, &quadtree.root, 0.5, params.softening, g);
+        for (i, &entity) in entities.iter().enumerate() {
+            let bh_force = calculate_barnes_hut_force(
+                entity,
+                bodies[i].0,
+                bodies[i].1,
+                &quadtree.root,
+                0.5,
+                params.softening,
+                g,
+            );
             let diff = (bh_force - brute_forces[i]).length();
             assert!(diff < 1.0, "BH force differs from brute force by {}", diff);
         }
diff --git a/crates/forces/src/core/mod.rs b/crates/forces/src/core/mod.rs
index 2da30ec..450f19c 100644
--- a/crates/forces/src/core/mod.rs
+++ b/crates/forces/src/core/mod.rs
@@ -15,9 +15,12 @@ pub mod prelude {
 
     // Re-export from newton_laws module
     pub use crate::core::newton_laws::{
-        AppliedForce, Distance, ForceImpulse, ForcesDiagnostics, ForcesDiagnosticsPlugin,
-        IntegratorKind, Mass, NewtonLawsPlugin, Norm, PairedForce, PairedForceInteraction,
-        Velocity, calculate_kinetic_energy, calculate_momentum, integrate_newton_second_law,
-        integrate_positions_symplectic_euler, update_forces_diagnostics,
+        AppliedForce, AppliedTorque, Distance, ForceImpulse, ForcesDiagnostics,
+        ForcesDiagnosticsPlugin, IntegratorKind, Mass, MomentOfInertia, NewtonLawsPlugin, Norm,
+        PairedForce, PairedForceInteraction, RotationalWorkEvent, Velocity, WorkDoneEvent,
+        calculate_angular_momentum, calculate_kinetic_energy, calculate_momentum,
+        calculate_rotational_kinetic_energy, calculate_torque_from_force,
+        integrate_newton_second_law, integrate_positions_symplectic_euler, integrate_torques,
+        update_forces_diagnostics,
     };
 }
diff --git a/crates/forces/src/core/newton_laws.rs b/crates/forces/src/core/newton_laws.rs
index 6188500..04e6e3f 100644
--- a/crates/forces/src/core/newton_laws.rs
+++ b/crates/forces/src/core/newton_laws.rs
@@ -1,3 +1,34 @@
+//! Entity-based Newtonian mechanics for non-continuum objects.
+//!
+//! # Architecture: Dual Physics Backend
+//!
+//! LP uses a **dual physics backend** to handle different types of matter:
+//!
+//! ## Entity-Based Physics (This Module)
+//! - **Scope**: Non-continuum objects (cameras, UI, large rigid bodies)
+//! - **Components**: `Mass`, `Velocity`, `AppliedForce`, `MomentOfInertia`, `AppliedTorque`
+//! - **Systems**: `integrate_newton_second_law`, `integrate_torques`
+//! - **Use cases**: Objects that don't deform or need particle-level resolution
+//!
+//! ## MPM Physics (Future: `crates/systems/mpm/`)
+//! - **Scope**: All continuum matter (water, soil, flesh, deformable solids)
+//! - **Representation**: Particles + grid (Material Point Method)
+//! - **Integration**: P2G transfer → grid solve → G2P transfer
+//! - **Use cases**: Anything that flows, deforms, or needs material-level physics
+//!
+//! ## Unified Conservation
+//! Both backends feed the same **energy ledger** (`crates/energy/conservation.rs`):
+//! - Entity backend: `WorkDoneEvent`, `RotationalWorkEvent` → ledger
+//! - MPM backend: Grid work events → same ledger (when implemented)
+//! - Diagnostics: Aggregate across both backends for global conservation tracking
+//!
+//! ## TODO (MPM Implementation)
+//! When MPM is implemented:
+//! - [ ] MPM will compute `L = Σ(r × m·v)` from particles (not `MomentOfInertia` components)
+//! - [ ] MPM will emit work events from grid-particle transfers
+//! - [ ] `ForcesDiagnostics` will aggregate both entity and MPM contributions
+//! - [ ] Gravity and other forces will have MPM-specific implementations
+
 use crate::PhysicsSet;
 use bevy::prelude::*;
 
@@ -175,8 +206,129 @@ pub struct Velocity {
     pub angvel: Vec3,
 }
 
+/// Component for moment of inertia (rotational mass).
+///
+/// **PHYSICS**: Moment of inertia I determines resistance to rotational acceleration.
+/// - For point mass: I = m·r² (kg·m²)
+/// - For rigid body: I = ∫ r²·dm (depends on mass distribution)
+///
+/// **UNITS**: kg·m² (kilogram-meter-squared)
+///
+/// **NOTE**: This is scalar moment of inertia for 2D rotation or principal axis in 3D.
+/// Full 3D rigid body needs inertia tensor (3×3 matrix), but most LP use cases are 2D or symmetric.
+#[derive(Component, Debug, Clone, Copy, Reflect)]
+#[reflect(Component)]
+pub struct MomentOfInertia {
+    /// Moment of inertia in kg·m²
+    pub value: f32,
+    /// Whether this object has infinite inertia (cannot rotate)
+    pub is_infinite: bool,
+}
+
+impl MomentOfInertia {
+    pub fn new(value: f32) -> Self {
+        debug_assert!(value > 0.0, "Moment of inertia cannot be negative or zero");
+        Self {
+            value: value.max(0.0001), // Prevent zero inertia
+            is_infinite: false,
+        }
+    }
+
+    pub fn infinite() -> Self {
+        Self {
+            value: f32::MAX,
+            is_infinite: true,
+        }
+    }
+
+    pub fn inverse(&self) -> f32 {
+        if self.is_infinite {
+            0.0
+        } else {
+            1.0 / self.value.max(f32::EPSILON)
+        }
+    }
+
+    /// Moment of inertia for a uniform disk/cylinder rotating about its axis
+    /// I = 0.5 * m * r²
+    pub fn disk(mass: f32, radius: f32) -> Self {
+        Self::new(0.5 * mass * radius * radius)
+    }
+
+    /// Moment of inertia for a solid sphere rotating about its center
+    /// I = 0.4 * m * r²
+    pub fn sphere(mass: f32, radius: f32) -> Self {
+        Self::new(0.4 * mass * radius * radius)
+    }
+
+    /// Moment of inertia for a point mass at distance r from axis
+    /// I = m * r²
+    pub fn point_mass(mass: f32, radius: f32) -> Self {
+        Self::new(mass * radius * radius)
+    }
+}
+
+/// Component representing torque applied to an entity.
+///
+/// **PHYSICS**: Torque τ = r × F (cross product of lever arm and force)
+/// - τ: Torque (N·m)
+/// - r: Position vector from rotation axis to force application point (m)
+/// - F: Applied force (N)
+///
+/// **UNITS**: Newton-meters (N·m) = kg·m²/s²
+///
+/// **ROTATIONAL ANALOG**: Torque is to angular acceleration what force is to linear acceleration.
+/// τ = I·α (where α is angular acceleration rad/s²)
+#[derive(Component, Debug, Clone, Reflect)]
+#[reflect(Component)]
+pub struct AppliedTorque {
+    /// Torque vector in Newton-meters (right-hand rule: thumb = axis, fingers = rotation direction)
+    pub torque: Vec3,
+    /// Duration the torque is applied (None for continuous)
+    pub duration: Option<f32>,
+    /// Elapsed time since torque began
+    pub elapsed: f32,
+}
+
+impl Default for AppliedTorque {
+    fn default() -> Self {
+        Self::new(Vec3::ZERO)
+    }
+}
+
+impl AppliedTorque {
+    pub fn new(torque: Vec3) -> Self {
+        Self {
+            torque,
+            duration: None,
+            elapsed: 0.0,
+        }
+    }
+
+    pub fn with_duration(mut self, duration: f32) -> Self {
+        self.duration = Some(duration);
+        self
+    }
+
+    pub fn is_expired(&self) -> bool {
+        if let Some(duration) = self.duration {
+            self.elapsed >= duration
+        } else {
+            false
+        }
+    }
+}
+
 /// Integrate Newton's Second Law to update velocities from forces.
 ///
+/// **ARCHITECTURE NOTE**: This is the **entity-based physics backend** for non-continuum objects
+/// (cameras, UI, large non-deformable bodies). For continuum matter (water, soil, deformables),
+/// LP uses MPM (Material Point Method) with particle-based integration.
+///
+/// **TODO (MPM)**: When MPM is implemented, continuum matter will use grid-based force integration
+/// (P2G transfer, grid solve, G2P transfer) instead of this entity-based approach. Both backends
+/// will emit work events to the same unified energy ledger.
+///
 /// **PHYSICS**: Newton's Second Law F = ma → a = F/m
 /// - F: Net force (Newtons)
 /// - m: Mass (kilograms)
@@ -194,18 +346,19 @@ pub struct Velocity {
 /// - Acceleration clamp: max 1000 m/s² (numerical stability, not IRL physics)
 ///
 /// **CONSERVATION**: Momentum conserved if all forces are action-reaction pairs (F_ab = -F_ba).
-/// Energy NOT conserved in this system (work tracking missing, no PE accounting).
+/// Energy conservation: work done by forces reported via WorkDoneEvent for ledger tracking.
 ///
 /// **NUMERICAL STABILITY**: Clamps extremely high accelerations to prevent explosive forces
 /// from breaking integration. If forces exceed limit, likely a bug (check softening, timestep, scales).
 pub fn integrate_newton_second_law(
     time: Res<Time>,
     config: Res<IntegrationConfig>,
-    mut query: Query<(&Mass, &mut Velocity, &mut AppliedForce)>,
+    mut query: Query<(Entity, &Mass, &mut Velocity, &mut AppliedForce)>,
+    mut work_events: MessageWriter<WorkDoneEvent>,
 ) {
     let dt = time.delta_secs();
 
-    for (mass, mut velocity, mut force) in query.iter_mut() {
+    for (entity, mass, mut velocity, mut force) in query.iter_mut() {
         if mass.is_infinite || mass.is_negligible() {
             continue;
         }
@@ -230,18 +383,28 @@ pub fn integrate_newton_second_law(
             acceleration
         };
 
+        // Calculate work using average velocity for better accuracy
+        // (Landau Vol 1 §6: dT = F·dr, where dr = v_avg·dt)
+        let v_old = velocity.linvel;
         velocity.linvel += acceleration * dt;
+        let v_avg = (v_old + velocity.linvel) * 0.5;
+        let work_done = force.force.dot(v_avg) * dt;
+
         force.elapsed += dt;
 
+        // Report work done for energy conservation tracking
+        if work_done.abs() > f32::EPSILON {
+            work_events.write(WorkDoneEvent {
+                entity,
+                work: work_done,
+            });
+        }
+
         // Clear accumulated force so subsequent systems can rebuild it per-frame
         force.force = Vec3::ZERO;
     }
 }
 
-// TODO: Work/energy reconciliation.
-// Work done by applied forces (∑ F·Δx) is not reported to the energy ledger yet.
-// When energy accounting is wired, emit an EnergyTransferEvent to keep conservation consistent.
-
 /// System to apply symplectic Euler integration for position updates
 pub fn integrate_positions_symplectic_euler(
     time: Res<Time>,
@@ -277,16 +440,147 @@ fn use_symplectic_euler(kind: Res<IntegratorKind>) -> bool {
     *kind == IntegratorKind::SymplecticEuler
 }
 
+/// Integrate torques to update angular velocities.
+///
+/// **ARCHITECTURE NOTE**: This is the **entity-based rigid body rotation** backend for non-continuum
+/// objects with discrete rotation (cameras, UI, large non-deformable bodies). For continuum matter,
+/// LP uses MPM where angular momentum emerges from particle positions: L = Σ(r_i × m_i·v_i).
+///
+/// **TODO (MPM)**: When MPM is implemented, continuum matter will compute angular momentum from
+/// particle distributions rather than rigid body components. MPM will emit rotational work events
+/// from grid-particle interactions to the same unified energy ledger.
+///
+/// **PHYSICS**: Rotational analog of Newton's Second Law: τ = I·α → α = τ/I
+/// - τ: Torque (N·m)
+/// - I: Moment of inertia (kg·m²)
+/// - α: Angular acceleration (rad/s²)
+/// - Δω = α·Δt (angular velocity update via Euler integration)
+///
+/// **UNITS**:
+/// - Torque: Newton-meters (N·m) = kg·m²/s²
+/// - Moment of inertia: kg·m²
+/// - Angular velocity: radians/second (rad/s)
+/// - Time: seconds (s)
+///
+/// **APPROXIMATIONS**:
+/// - Explicit Euler: ω_new = ω_old + α·dt (first-order accurate)
+/// - Assumes principal axis rotation (scalar moment of inertia)
+/// - Full 3D requires inertia tensor (future work)
+///
+/// **CONSERVATION**: Angular momentum conserved if all torques are balanced.
+/// Rotational work reported via RotationalWorkEvent for ledger tracking.
+pub fn integrate_torques(
+    time: Res<Time>,
+    config: Res<IntegrationConfig>,
+    mut query: Query<(Entity, &MomentOfInertia, &mut Velocity, &mut AppliedTorque)>,
+    mut rotational_work_events: MessageWriter<RotationalWorkEvent>,
+) {
+    let dt = time.delta_secs();
+
+    for (entity, inertia, mut velocity, mut torque) in query.iter_mut() {
+        if inertia.is_infinite {
+            continue; // Cannot rotate
+        }
+
+        if torque.is_expired() {
+            torque.torque = Vec3::ZERO;
+            continue;
+        }
+
+        // Calculate angular acceleration: α = τ/I
+        let angular_acceleration = torque.torque * inertia.inverse();
+
+        if !angular_acceleration.is_finite() {
+            torque.torque = Vec3::ZERO;
+            continue;
+        }
+
+        // Cap extremely high angular accelerations for stability
+        // Using same threshold as linear (1000 rad/s² ≈ 57000°/s²)
+        let max_angular_acceleration = config.max_acceleration; // rad/s²
+        let angular_acceleration = if angular_acceleration.norm_squared()
+            > max_angular_acceleration * max_angular_acceleration
+        {
+            angular_acceleration.normalize() * max_angular_acceleration
+        } else {
+            angular_acceleration
+        };
+
+        // Calculate rotational work using average angular velocity
+        // W_rot = τ·Δθ where Δθ = ω_avg·dt (angle rotated)
+        let omega_old = velocity.angvel;
+        velocity.angvel += angular_acceleration * dt;
+        let omega_avg = (omega_old + velocity.angvel) * 0.5;
+        let delta_theta = omega_avg * dt; // Angle rotated (radians)
+        let work_done = torque.torque.dot(delta_theta);
+
+        torque.elapsed += dt;
+
+        // Report rotational work done for energy conservation tracking
+        if work_done.abs() > f32::EPSILON {
+            rotational_work_events.write(RotationalWorkEvent {
+                entity,
+                work: work_done,
+            });
+        }
+
+        // Clear accumulated torque
+        torque.torque = Vec3::ZERO;
+    }
+}
+
 /// Calculate momentum of an object
 pub fn calculate_momentum(mass: &Mass, velocity: &Velocity) -> Vec3 {
     mass.value * velocity.linvel
 }
 
-/// Calculate kinetic energy of an object
+/// Calculate kinetic energy of an object (translational only)
 pub fn calculate_kinetic_energy(mass: &Mass, velocity: &Velocity) -> f32 {
     0.5 * mass.value * velocity.linvel.norm_squared()
 }
 
+/// Calculate angular momentum of an object.
+///
+/// **PHYSICS**: L = I·ω (angular momentum = moment of inertia × angular velocity)
+/// - L: Angular momentum (kg·m²/s)
+/// - I: Moment of inertia (kg·m²)
+/// - ω: Angular velocity (rad/s)
+///
+/// **UNITS**: kg·m²/s
+///
+/// **CONSERVATION**: L is conserved when no external torques act on the system.
+pub fn calculate_angular_momentum(inertia: &MomentOfInertia, velocity: &Velocity) -> Vec3 {
+    inertia.value * velocity.angvel
+}
+
+/// Calculate rotational kinetic energy of an object.
+///
+/// **PHYSICS**: KE_rot = 0.5·I·ω² (rotational kinetic energy)
+/// - KE_rot: Rotational kinetic energy (J)
+/// - I: Moment of inertia (kg·m²)
+/// - ω: Angular velocity (rad/s)
+///
+/// **UNITS**: Joules (J) = kg·m²/s²
+///
+/// **TOTAL KINETIC ENERGY**: KE_total = KE_trans + KE_rot = 0.5·m·v² + 0.5·I·ω²
+pub fn calculate_rotational_kinetic_energy(inertia: &MomentOfInertia, velocity: &Velocity) -> f32 {
+    0.5 * inertia.value * velocity.angvel.norm_squared()
+}
+
+/// Calculate torque from an off-center force.
+///
+/// **PHYSICS**: τ = r × F (cross product of lever arm and force)
+/// - τ: Torque (N·m)
+/// - r: Position vector from rotation axis to force application point (m)
+/// - F: Force vector (N)
+///
+/// **RIGHT-HAND RULE**: Fingers point from r to F, thumb points along τ
+///
+/// **UNITS**: Newton-meters (N·m) = kg·m²/s²
+pub fn calculate_torque_from_force(lever_arm: Vec3, force: Vec3) -> Vec3 {
+    lever_arm.cross(force)
+}
+
 /// Represents a pair of entities for force calculations (Newton's Third Law)
 #[derive(Debug, Clone, Copy)]
 pub struct ForcePair<'a> {
@@ -324,6 +618,39 @@ impl ForceImpulse {
     }
 }
 
+/// Event for reporting work done by forces (W = F·dx)
+///
+/// **ARCHITECTURE**: Part of unified conservation interface. Entity-based physics emits these
+/// events when forces do work on entities. Energy crate (`conservation.rs`) listens and records
+/// to the energy ledger.
+///
+/// **TODO (MPM)**: When MPM is implemented, it will emit equivalent work events from grid-particle
+/// transfers. Both entity and MPM events feed the same `EnergyBalance` ledger for unified tracking.
+///
+/// Energy crate listens to this to track kinetic energy changes.
+#[derive(Message)]
+pub struct WorkDoneEvent {
+    pub entity: Entity,
+    pub work: f32, // Joules
+}
+
+/// Event for reporting rotational work done by torques (W_rot = τ·Δθ)
+///
+/// **ARCHITECTURE**: Part of unified conservation interface. Entity-based physics emits these
+/// events when torques do rotational work on rigid bodies. Energy crate (`conservation.rs`)
+/// listens and records to the energy ledger.
+///
+/// **TODO (MPM)**: When MPM is implemented, rotational work will emerge from particle motions
+/// around centers of mass. MPM will emit equivalent events computed from particle angular momentum
+/// changes. Both entity and MPM events feed the same `EnergyBalance` ledger.
+///
+/// Energy crate listens to this to track rotational kinetic energy changes.
+#[derive(Message)]
+pub struct RotationalWorkEvent {
+    pub entity: Entity,
+    pub work: f32, // Joules
+}
+
 /// Plugin that adds Newton's Laws mechanics systems in the correct order
 #[derive(Default)]
 pub struct NewtonLawsPlugin;
@@ -335,6 +662,8 @@ impl Plugin for NewtonLawsPlugin {
         app.init_resource::<IntegratorKind>()
             .init_resource::<IntegrationConfig>()
             .add_message::<ForceImpulse>()
+            .add_message::<WorkDoneEvent>()
+            .add_message::<RotationalWorkEvent>()
             // Configure physics sets
             .configure_sets(
                 Update,
@@ -345,10 +674,14 @@ impl Plugin for NewtonLawsPlugin {
                 )
                     .chain(),
             )
-            // Apply forces and integrate - chain integrate_newton_second_law and apply_impulses
+            // Apply forces and torques, then integrate - chain all three
             .add_systems(
                 Update,
-                (integrate_newton_second_law, apply_impulses)
+                (
+                    integrate_newton_second_law,
+                    integrate_torques,
+                    apply_impulses,
+                )
                     .chain()
                     .in_set(PhysicsSet::ApplyForces),
             )
@@ -404,32 +737,48 @@ pub fn apply_impulses(
     }
 }
 
-/// Snapshot diagnostics for total momentum and kinetic energy.
+/// Snapshot diagnostics for total momentum, angular momentum, and kinetic energy.
 #[derive(Resource, Debug, Default, Clone)]
 pub struct ForcesDiagnostics {
     pub total_momentum: Vec3,
+    pub total_angular_momentum: Vec3,
     pub total_kinetic_energy: f32,
+    pub total_rotational_kinetic_energy: f32,
 }
 
 /// Updates diagnostics after velocity changes are applied.
 pub fn update_forces_diagnostics(
     mut diagnostics: ResMut<ForcesDiagnostics>,
-    query: Query<(&Mass, &Velocity)>,
+    query: Query<(&Mass, &Velocity, Option<&MomentOfInertia>)>,
 ) {
     let mut total_momentum = Vec3::ZERO;
+    let mut total_angular_momentum = Vec3::ZERO;
     let mut total_kinetic_energy = 0.0;
+    let mut total_rotational_kinetic_energy = 0.0;
 
-    for (mass, velocity) in &query {
+    for (mass, velocity, inertia_opt) in &query {
         if mass.is_infinite {
             continue;
         }
 
+        // Linear momentum and kinetic energy
         total_momentum += calculate_momentum(mass, velocity);
         total_kinetic_energy += calculate_kinetic_energy(mass, velocity);
+
+        // Angular momentum and rotational kinetic energy (if entity has moment of inertia)
+        if let Some(inertia) = inertia_opt {
+            if !inertia.is_infinite {
+                total_angular_momentum += calculate_angular_momentum(inertia, velocity);
+                total_rotational_kinetic_energy +=
+                    calculate_rotational_kinetic_energy(inertia, velocity);
+            }
+        }
     }
 
     diagnostics.total_momentum = total_momentum;
+    diagnostics.total_angular_momentum = total_angular_momentum;
     diagnostics.total_kinetic_energy = total_kinetic_energy;
+    diagnostics.total_rotational_kinetic_energy = total_rotational_kinetic_energy;
 }
 
 /// Plugin to enable Newton-law diagnostics.
@@ -509,4 +858,80 @@ mod tests {
         // KE = 0.5 * m * v² = 0.5 * 2.0 * (3² + 4²) = 0.5 * 2.0 * 25 = 25.0
         assert_eq!(ke, 25.0);
     }
+
+    #[test]
+    fn test_angular_momentum_calculation() {
+        // L = I·ω (angular momentum = moment of inertia × angular velocity)
+        let inertia = MomentOfInertia::new(2.0); // kg·m²
+        let velocity = Velocity {
+            linvel: Vec3::ZERO,
+            angvel: Vec3::new(0.0, 0.0, 5.0), // 5 rad/s around Z axis
+        };
+        let angular_momentum = calculate_angular_momentum(&inertia, &velocity);
+        // L = I·ω = 2.0 × 5.0 = 10.0 kg·m²/s
+        assert_eq!(angular_momentum, Vec3::new(0.0, 0.0, 10.0));
+    }
+
+    #[test]
+    fn test_rotational_kinetic_energy_calculation() {
+        // KE_rot = 0.5·I·ω²
+        let inertia = MomentOfInertia::new(2.0); // kg·m²
+        let velocity = Velocity {
+            linvel: Vec3::ZERO,
+            angvel: Vec3::new(0.0, 0.0, 3.0), // 3 rad/s around Z axis
+        };
+        let ke_rot = calculate_rotational_kinetic_energy(&inertia, &velocity);
+        // KE_rot = 0.5 × 2.0 × 3² = 0.5 × 2.0 × 9 = 9.0 J
+        assert_eq!(ke_rot, 9.0);
+    }
+
+    #[test]
+    fn test_torque_from_force_calculation() {
+        // τ = r × F (cross product)
+        // If force is applied at distance r from rotation axis, creates torque
+        let lever_arm = Vec3::new(2.0, 0.0, 0.0); // 2m along X axis
+        let force = Vec3::new(0.0, 5.0, 0.0); // 5N along Y axis
+        let torque = calculate_torque_from_force(lever_arm, force);
+        // τ = r × F = (2,0,0) × (0,5,0) = (0,0,10) N·m
+        // Right-hand rule: fingers point from X to Y, thumb points along Z
+        assert_eq!(torque, Vec3::new(0.0, 0.0, 10.0));
+    }
+
+    #[test]
+    fn test_moment_of_inertia_formulas() {
+        let mass = 10.0; // kg
+        let radius = 2.0; // m
+
+        // Disk: I = 0.5·m·r²
+        let disk = MomentOfInertia::disk(mass, radius);
+        assert_eq!(disk.value, 0.5 * 10.0 * 2.0 * 2.0); // = 20.0
+
+        // Sphere: I = 0.4·m·r² (actually 2/5)
+        let sphere = MomentOfInertia::sphere(mass, radius);
+        assert_eq!(sphere.value, 0.4 * 10.0 * 2.0 * 2.0); // = 16.0
+
+        // Point mass: I = m·r²
+        let point = MomentOfInertia::point_mass(mass, radius);
+        assert_eq!(point.value, 10.0 * 2.0 * 2.0); // = 40.0
+    }
+
+    #[test]
+    fn test_total_energy_includes_rotation() {
+        // Total KE = translational + rotational
+        let mass = Mass::new(2.0);
+        let inertia = MomentOfInertia::new(3.0);
+        let velocity = Velocity {
+            linvel: Vec3::new(4.0, 0.0, 0.0), // 4 m/s
+            angvel: Vec3::new(0.0, 0.0, 2.0), // 2 rad/s
+        };
+
+        let ke_trans = calculate_kinetic_energy(&mass, &velocity);
+        // KE_trans = 0.5 × 2.0 × 4² = 16.0 J
+
+        let ke_rot = calculate_rotational_kinetic_energy(&inertia, &velocity);
+        // KE_rot = 0.5 × 3.0 × 2² = 6.0 J
+
+        let total_ke = ke_trans + ke_rot;
+        assert_eq!(total_ke, 22.0); // 16.0 + 6.0 = 22.0 J
+    }
 }
diff --git a/crates/forces/src/lib.rs b/crates/forces/src/lib.rs
index 788c4fd..8dc3d3a 100644
--- a/crates/forces/src/lib.rs
+++ b/crates/forces/src/lib.rs
@@ -30,15 +30,13 @@ pub struct ForcesPlugin;
 
 impl Plugin for ForcesPlugin {
     fn build(&self, app: &mut App) {
-        app.add_plugins(NewtonLawsPlugin)
+        app.add_plugins((NewtonLawsPlugin, core::gravity::GravityPlugin::new()))
             .register_type::<core::newton_laws::Mass>()
             .register_type::<core::newton_laws::Velocity>()
             .register_type::<core::newton_laws::AppliedForce>()
             .register_type::<core::gravity::GravityAffected>()
             .register_type::<core::gravity::GravitySource>()
-            .register_type::<core::gravity::MassiveBody>()
-            .init_resource::<core::gravity::GravityParams>()
-            .init_resource::<core::gravity::UniformGravity>();
+            .register_type::<core::gravity::MassiveBody>();
     }
 }