diff --git a/docs/dc_plus/injection_outage.md b/docs/dc_plus/injection_outage.md
new file mode 100644
index 0000000..344ebbb
--- /dev/null
+++ b/docs/dc_plus/injection_outage.md
@@ -0,0 +1,32 @@
+# Injection Outage
+
+`non_voltage_regulating_injection_outage_dx` computes one-step AC state updates for a batch of injection outages around the N-0 operating point.
+
+The Jacobian is kept fixed. Each outage only contributes an active and reactive power mismatch at the affected bus, and the batched update is evaluated with the base inverse Jacobian.
+
+$$
+\Delta X = - M J^{-T}
+$$
+
+where `J` is the base-case Jacobian and `M` is the batch of outage mismatch vectors in Jacobian ordering.
+
+## Batch API
+
+The function expects a 2D outage array with shape `(n_contingencies, n_outages)`.
+
+- Each row is one contingency.
+- A row may contain multiple outage indices.
+- Negative entries are ignored and can be used as padding.
+- The return value has shape `(n_contingencies, n_eq)`.
+
+## Scope
+
+This update is intended for non-voltage-regulating injections such as loads and other fixed-power injections.
+
+## Generator Behavior
+
+The function can also be applied to voltage-regulating generator outages, but this is only an approximation. A voltage-regulating generator outage should normally trigger a PV to PQ bus-type change. The current implementation does not model that switch, so it does not capture the resulting voltage drop at the regulated bus.
+
+## HVDC
+
+- HVDC converter outages should be outaged in pairs. DC+ models the HVDC line as a pair of generator + load injections, one at each terminal.
diff --git a/src/dc_plus/example_grids/pypowsbl/example_grids.py b/src/dc_plus/example_grids/pypowsbl/example_grids.py
index 2aae34b..5cbd517 100644
--- a/src/dc_plus/example_grids/pypowsbl/example_grids.py
+++ b/src/dc_plus/example_grids/pypowsbl/example_grids.py
@@ -983,6 +983,15 @@ def basic_node_breaker_network_powsybl() -> pypowsybl.network.Network:
     pypowsybl.network.create_metrix_tutorial_six_buses_network,  # HVDC
 ]
 
+POWSYBL_NETWORKS_SHORT_LIST = [
+    pypowsybl.network.create_ieee9,
+    pypowsybl.network.create_ieee14,
+    pypowsybl.network.create_ieee57,
+    pypowsybl.network.create_micro_grid_nl_network,
+    basic_node_breaker_network_powsybl,
+    create_complex_grid_battery_hvdc_svc_3w_trafo,
+]
+
 # TODO: check this network and add to POWSYBL_NETWORKS if possible (currently fails)
 POWSYBL_NETWORKS_NOT_IMPLEMENTED = [
     pypowsybl.network.create_micro_grid_be_network,
diff --git a/src/dc_plus/importing/import_schema.py b/src/dc_plus/importing/import_schema.py
index 5a4db84..81d10e0 100644
--- a/src/dc_plus/importing/import_schema.py
+++ b/src/dc_plus/importing/import_schema.py
@@ -114,6 +114,10 @@ class BusParamSchema(pa.DataFrameModel):
     voltage_magnitude: pat.Series[float] = pa.Field(coerce=True, nullable=True)
     voltage_angle: pat.Series[float] = pa.Field(coerce=True, nullable=True)
     bus_type: pat.Series[int] = pa.Field(coerce=True, description="0:Slack, 1:PV, 2:PQ")
+    is_angle_reference: pat.Series[bool] = pa.Field(
+        coerce=True,
+        description="True for the bus used as the angle reference, independent of its voltage-control mode.",
+    )
     grid_island_id: pat.Series[int] = pa.Field(
         coerce=True, description="ID of the grid island the bus belongs to. 0 indicates the main grid."
     )
diff --git a/src/dc_plus/importing/pandapower/pandapower_import.py b/src/dc_plus/importing/pandapower/pandapower_import.py
index 3337a70..ddf2c7c 100644
--- a/src/dc_plus/importing/pandapower/pandapower_import.py
+++ b/src/dc_plus/importing/pandapower/pandapower_import.py
@@ -272,11 +272,10 @@ def _get_buses_pandapower(net: pandapowerNet, slack_id: NotImplementedError) ->
             continue
 
         res_bus = net.res_bus.loc[idx]
+        is_angle_reference = idx == slack_id
 
         # Determine bus type
-        if idx == slack_id:
-            bus_type = BusType.SLACK
-        elif idx in net.gen["bus"].values or idx in net.sgen["bus"].values:
+        if idx in net.gen["bus"].values or idx in net.sgen["bus"].values:
             # Check if voltage control is enabled
             is_voltage_controlled = False
             if idx in net.gen["bus"].values:
@@ -285,8 +284,13 @@ def _get_buses_pandapower(net: pandapowerNet, slack_id: NotImplementedError) ->
             if not is_voltage_controlled and idx in net.sgen["bus"].values:
                 sgen_at_bus = net.sgen[net.sgen["bus"] == idx]
                 is_voltage_controlled = sgen_at_bus.get("controllable", pd.Series([False])).any()
+        else:
+            is_voltage_controlled = is_angle_reference
 
-            bus_type = BusType.PV if is_voltage_controlled else BusType.PQ
+        if is_angle_reference and is_voltage_controlled:
+            bus_type = BusType.SLACK
+        elif is_voltage_controlled:
+            bus_type = BusType.PV
         else:
             bus_type = BusType.PQ
 
@@ -302,6 +306,7 @@ def _get_buses_pandapower(net: pandapowerNet, slack_id: NotImplementedError) ->
                 "voltage_magnitude": res_bus["vm_pu"] if pd.notna(res_bus["vm_pu"]) else np.nan,
                 "voltage_angle": (np.deg2rad(res_bus["va_degree"]) if pd.notna(res_bus["va_degree"]) else np.nan),
                 "bus_type": int(bus_type),
+                "is_angle_reference": is_angle_reference,
                 "grid_island_id": grid_island_id,
             }
         )
diff --git a/src/dc_plus/importing/powsybl/powsybl_import.py b/src/dc_plus/importing/powsybl/powsybl_import.py
index 1d0b45e..7768b65 100644
--- a/src/dc_plus/importing/powsybl/powsybl_import.py
+++ b/src/dc_plus/importing/powsybl/powsybl_import.py
@@ -20,6 +20,7 @@
     LimitParamSchema,
     ShuntParamSchema,
 )
+from dc_plus.interfaces.network_information import BusType
 
 logger = logging.getLogger(__name__)
 
@@ -907,17 +908,27 @@ def _get_dangling_buses(net: Network) -> BusParamSchema:
     # TODO: get values depending on the connected state and the bus_id of the dangling line
     dangling_buses["grid_island_id"] = 0
     dangling_buses["bus_type"] = 2  # PQ bus
+    dangling_buses["is_angle_reference"] = False
     bus_order_int_ids = _get_bus_ids_with_dangling_buses(net)
     dangling_buses = dangling_buses.merge(bus_order_int_ids, how="left", on=["id_str"])
 
     dangling_buses = dangling_buses[
-        ["id_str", "id_int", "name", "voltage_magnitude", "voltage_angle", "bus_type", "grid_island_id"]
+        [
+            "id_str",
+            "id_int",
+            "name",
+            "voltage_magnitude",
+            "voltage_angle",
+            "bus_type",
+            "is_angle_reference",
+            "grid_island_id",
+        ]
     ]
     dangling_buses = BusParamSchema.validate(dangling_buses)
     return dangling_buses
 
 
-def _get_buses_powsybl(net: Network, slack_id: str, injections: InjectionParamSchema) -> BusParamSchema:
+def _get_buses_powsybl(net: Network, slack_id: str, injections: InjectionParamSchema, reference_id: str) -> BusParamSchema:
     """Get the bus parameters of the network.
 
     Gets the bus parameters from the network.
@@ -930,6 +941,8 @@ def _get_buses_powsybl(net: Network, slack_id: str, injections: InjectionParamSc
         The id of the slack bus.
     injections : InjectionParamSchema
         The injections of the network.
+    reference_id : str
+        The id of the reference bus for angle reference.
 
     Returns
     -------
@@ -942,6 +955,7 @@ def _get_buses_powsybl(net: Network, slack_id: str, injections: InjectionParamSc
     buses.reset_index(drop=False, inplace=True)
     buses.rename(columns={"id": "id_str"}, inplace=True)
     buses["bus_type"] = 2  # PQ bus is default
+    buses["is_angle_reference"] = buses["id_str"].eq(reference_id)
 
     buses.rename(columns={"v_mag": "voltage_magnitude", "v_angle": "voltage_angle"}, inplace=True)
 
@@ -966,15 +980,29 @@ def _get_buses_powsybl(net: Network, slack_id: str, injections: InjectionParamSc
     buses = buses.merge(bus_order_int_ids, how="left", on=["id_str"])
 
     # now the main grid "0" has only connected_component == 0 and synchronous_component == 0
-    buses = buses[["id_str", "id_int", "name", "voltage_magnitude", "voltage_angle", "bus_type", "grid_island_id"]]
+    buses = buses[
+        [
+            "id_str",
+            "id_int",
+            "name",
+            "voltage_magnitude",
+            "voltage_angle",
+            "bus_type",
+            "is_angle_reference",
+            "grid_island_id",
+        ]
+    ]
 
     dangling_buses = _get_dangling_buses(net)
     buses = pd.concat([buses, dangling_buses], ignore_index=True)
     buses.sort_values(by=["id_int"], inplace=True)
-    # set bus types
-    pv_buses = injections[(injections["voltage_regulation"])]["regulated_bus_id_int"].unique()
-    buses.loc[buses["id_int"].isin(pv_buses), "bus_type"] = 1  # PV bus
-    buses.loc[buses["id_str"] == slack_id, "bus_type"] = 0  # slack bus
+    # Bus type tracks voltage control; the reference-angle role is stored separately.
+    voltage_control_buses = injections[(injections["connected"]) & (injections["voltage_regulation"])][
+        "regulated_bus_id_int"
+    ].unique()
+    buses.loc[buses["id_int"].isin(voltage_control_buses), "bus_type"] = BusType.PV
+    slack_mask = buses["id_str"].eq(slack_id) & buses["id_int"].isin(voltage_control_buses)
+    buses.loc[slack_mask, "bus_type"] = BusType.SLACK
 
     buses = BusParamSchema.validate(buses)
     return buses
diff --git a/src/dc_plus/importing/powsybl/powsybl_import_helpers.py b/src/dc_plus/importing/powsybl/powsybl_import_helpers.py
index 8e166bf..e3498b0 100644
--- a/src/dc_plus/importing/powsybl/powsybl_import_helpers.py
+++ b/src/dc_plus/importing/powsybl/powsybl_import_helpers.py
@@ -9,7 +9,7 @@
 
 import pypowsybl
 import pypowsybl.loadflow
-from pypowsybl.loadflow import ConnectedComponentMode, VoltageInitMode
+from pypowsybl.loadflow import ComponentMode, VoltageInitMode
 from pypowsybl.loadflow import Parameters as LoadFlowParameters
 from pypowsybl.network import Network
 
@@ -56,7 +56,7 @@ def select_a_generator_as_slack_and_run_loadflow(network: Network) -> None:
         read_slack_bus=True,
         distributed_slack=True,
         use_reactive_limits=True,
-        connected_component_mode=ConnectedComponentMode.MAIN,  # ConnectedComponentMode
+        component_mode=ComponentMode.MAIN_CONNECTED,
     )
 
     loadflow_res = pypowsybl.loadflow.run_ac(network, parameters=powsybl_loadflow_parameter)[0]
diff --git a/src/dc_plus/importing/powsybl/powsybl_loadflow_parameter.py b/src/dc_plus/importing/powsybl/powsybl_loadflow_parameter.py
index f447bff..e2bc69f 100644
--- a/src/dc_plus/importing/powsybl/powsybl_loadflow_parameter.py
+++ b/src/dc_plus/importing/powsybl/powsybl_loadflow_parameter.py
@@ -10,7 +10,7 @@
 import logging
 from typing import Literal
 
-from pypowsybl.loadflow import ConnectedComponentMode, VoltageInitMode
+from pypowsybl.loadflow import ComponentMode, VoltageInitMode
 from pypowsybl.loadflow import Parameters as LoadFlowParameters
 
 logger = logging.getLogger(__name__)
@@ -59,7 +59,7 @@ def get_powsybl_loadflow_parameter(
         # balance_type=BalanceType.PROPORTIONAL_TO_GENERATION_REMAINING_MARGIN,  # BalanceType
         dc_use_transformer_ratio=True,
         countries_to_balance=None,  # Sequence[str]
-        connected_component_mode=ConnectedComponentMode.MAIN,  # ConnectedComponentMode
+        component_mode=ComponentMode.MAIN_CONNECTED,
         dc_power_factor=None,
         provider_parameters=provider_param,  # Dict[str, str]
     )
diff --git a/src/dc_plus/interfaces/jacobian_network_data.py b/src/dc_plus/interfaces/jacobian_network_data.py
index db6faf5..c64fd19 100644
--- a/src/dc_plus/interfaces/jacobian_network_data.py
+++ b/src/dc_plus/interfaces/jacobian_network_data.py
@@ -51,7 +51,7 @@ def _get_jacobian_data_from_network_data(dynamic_network_data: DynamicNetworkInf
         jacobian_index_in_use=jacobian_index_in_use,
         pointer_to_original_bus=np.arange(dynamic_network_data.n_buses, dtype=np.int32),
         jacobian=jacobian,
-        inverse_jacobian=sparse_inv(jacobian).toarray(),
+        inverse_jacobian=sparse_inv(jacobian.tocsc()).toarray(),
         is_angle_component=is_angle_component,
         is_magnitude_component=is_magnitude_component,
         pvpq_indices=pvpq_indices,
@@ -86,7 +86,7 @@ def get_jacobian_from_network_data(
     pqpv_indices = dynamic_network_data.pvpq_buses_indices_pvpq_order
 
     voltage = voltage_magnitudes * np.exp(1.0j * voltage_angles)
-    if np.any(voltage == 0.0):
+    if np.any(np.isclose(voltage, 0.0, atol=1e-10)):
         raise ValueError("Voltage magnitudes must be strictly positive to construct the Jacobian.")
 
     voltage_norm = voltage / abs(voltage)
diff --git a/src/dc_plus/interfaces/network_information.py b/src/dc_plus/interfaces/network_information.py
index 0301d4b..0db7d39 100644
--- a/src/dc_plus/interfaces/network_information.py
+++ b/src/dc_plus/interfaces/network_information.py
@@ -349,6 +349,14 @@ class DynamicNetworkInformation:
     2: pq
     """
 
+    bus_is_angle_reference: Bool[np.ndarray, " n_buses"]
+    """Indicates which bus is used as the angle reference.
+
+    This is distinct from the bus voltage-control mode. A reference bus can lose
+    voltage control after an outage and become a PQ bus while still fixing the
+    global angle reference.
+    """
+
     # injection data
 
     injection_to_bus: Int[np.ndarray, " n_injections n_timestep"]
@@ -402,14 +410,14 @@ def n_branches(self) -> int:
 
     @property
     def slack_indices(self) -> np.ndarray:
-        """Return the indices of the slack bus.
+        """Return the indices of the angle-reference bus.
 
         Returns
         -------
         slack_index : np.ndarray
-            The indices of the slack bus.
+            The indices of the angle-reference bus.
         """
-        slack_index = np.flatnonzero(self.bus_type == BusType.SLACK)
+        slack_index = np.flatnonzero(self.bus_is_angle_reference)
         return slack_index
 
     def is_pv_bus(self, bus_index: int) -> bool:
@@ -488,7 +496,7 @@ def pv_buses_mask(self) -> np.ndarray:
         pv_buses_mask : np.ndarray
             A boolean mask indicating which buses are PV buses.
         """
-        pv_buses_mask = self.bus_type == BusType.PV
+        pv_buses_mask = (self.bus_type == BusType.PV) & ~self.bus_is_angle_reference
         return pv_buses_mask
 
     @property
@@ -505,14 +513,14 @@ def pq_buses_mask(self) -> np.ndarray:
 
     @property
     def pvpq_buses_mask(self) -> np.ndarray:
-        """Return a boolean mask indicating which buses are PV or PQ buses.
+        """Return a boolean mask indicating which buses participate in angle equations.
 
         Returns
         -------
         pvpq_buses_mask : np.ndarray
             A boolean mask indicating which buses are PV or PQ buses.
         """
-        pvpq_buses_mask = (self.bus_type == BusType.PV) | (self.bus_type == BusType.PQ)
+        pvpq_buses_mask = ~self.bus_is_angle_reference
         return pvpq_buses_mask
 
     @property
@@ -524,7 +532,7 @@ def pv_buses_indices(self) -> np.ndarray:
         pv_buses_indices : np.ndarray
             The indices of the PV buses.
         """
-        pv_buses_indices = np.flatnonzero(self.bus_type == BusType.PV)
+        pv_buses_indices = np.flatnonzero(self.pv_buses_mask)
         return pv_buses_indices
 
     @property
@@ -536,31 +544,32 @@ def pq_buses_indices(self) -> np.ndarray:
         pq_buses_indices : np.ndarray
             The indices of the PQ buses.
         """
-        pq_buses_indices = np.flatnonzero(self.bus_type == BusType.PQ)
+        pq_buses_indices = np.flatnonzero(self.pq_buses_mask)
         return pq_buses_indices
 
     @property
     def pvpq_buses_indices(self) -> np.ndarray:
-        """Return the indices of the PV and PQ buses.
+        """Return the indices of buses participating in angle equations.
 
         Returns
         -------
         pvpq_buses_indices : np.ndarray
             The indices of the PV and PQ buses.
         """
-        pvpq_buses_indices = np.flatnonzero((self.bus_type == BusType.PV) | (self.bus_type == BusType.PQ))
+        pvpq_buses_indices = np.flatnonzero(self.pvpq_buses_mask)
         return pvpq_buses_indices
 
     @property
     def pvpq_buses_indices_pvpq_order(self) -> np.ndarray:
-        """Return the indices of the PV and PQ buses.
+        """Return the indices of buses participating in angle equations.
 
         Returns
         -------
         pvpq_buses_indices : np.ndarray
             The indices of the PV and PQ buses.
         """
-        pvpq_buses_indices = np.concatenate((self.pv_buses_indices, self.pq_buses_indices))
+        pq_angle_indices = self.pq_buses_indices[~self.bus_is_angle_reference[self.pq_buses_indices]]
+        pvpq_buses_indices = np.concatenate((self.pv_buses_indices, pq_angle_indices))
         return pvpq_buses_indices
 
     # injections
@@ -691,6 +700,9 @@ def _check_network_data_consistency(
     assert dynamic_network_data.bus_voltage_magnitudes.shape[0] == dynamic_network_data.bus_type.shape[0], (
         "Inconsistent number of buses between bus_voltage_magnitudes and bus_type."
     )
+    assert dynamic_network_data.bus_voltage_magnitudes.shape[0] == dynamic_network_data.bus_is_angle_reference.shape[0], (
+        "Inconsistent number of buses between bus_voltage_magnitudes and bus_is_angle_reference."
+    )
 
     # check injection data
     assert dynamic_network_data.injection_to_bus.shape[0] == dynamic_network_data.injection_active_power.shape[0], (
diff --git a/src/dc_plus/jax/bsdf.py b/src/dc_plus/jax/bsdf.py
index a84f597..12ba4b7 100644
--- a/src/dc_plus/jax/bsdf.py
+++ b/src/dc_plus/jax/bsdf.py
@@ -80,8 +80,9 @@ def _compute_bsdf_update_impl(
     angle_component_indices: Int[jnp.ndarray, " n_eq_jacobian"],
     magnitude_component_indices: Int[jnp.ndarray, " n_eq_jacobian"],
 ) -> Float[jnp.ndarray, " n_eq n_eq"]:
-    """Legacy implementation: materialize updated inverse transpose."""
-    dtype = jacobian_inv_transposed.dtype
+    """Materialize updated inverse transpose."""
+    # get dtype -> either float32 or float64 depending on input
+    dtype_input = v_mag_hat.dtype
     n_eq = jacobian_inv_transposed.shape[0]
 
     branch_from_old = jnp.take(branch_from, branches_connected_to_bus_b, axis=0)
@@ -202,7 +203,7 @@ def _gather_positions(indices: jnp.ndarray) -> Tuple[jnp.ndarray, jnp.ndarray]:
     )
 
     k_shape = (k_max, k_max)
-    delta_block = jnp.zeros(k_shape, dtype=dtype)
+    delta_block = jnp.zeros(k_shape, dtype=dtype_input)
 
     def _accumulate(
         delta: jnp.ndarray, pos: jnp.ndarray, valid: jnp.ndarray, weight: float, target: jnp.ndarray
@@ -236,8 +237,8 @@ def _accumulate(
     theta_row = jnp.where(theta_mask, theta_new_pos, 0)
     mag_row = jnp.where(mag_mask, mag_new_pos, 0)
 
-    minus_one = jnp.asarray(-1.0, dtype=dtype)
-    zero_val = jnp.asarray(0.0, dtype=dtype)
+    minus_one = jnp.asarray(-1.0, dtype=dtype_input)
+    zero_val = jnp.asarray(0.0, dtype=dtype_input)
     theta_update = jnp.where(theta_mask, minus_one, zero_val)
     mag_update = jnp.where(mag_mask, minus_one, zero_val)
     delta_block = delta_block.at[theta_row, theta_row].add(theta_update)
diff --git a/src/dc_plus/jax/injection_outage.py b/src/dc_plus/jax/injection_outage.py
new file mode 100644
index 0000000..d2ae4a6
--- /dev/null
+++ b/src/dc_plus/jax/injection_outage.py
@@ -0,0 +1,467 @@
+# Copyright 2026 50Hertz Transmission GmbH and Elia Transmission Belgium SA/NV
+#
+# This Source Code Form is subject to the terms of the Mozilla Public License, v. 2.0.
+# If a copy of the MPL was not distributed with this file,
+# you can obtain one at https://mozilla.org/MPL/2.0/.
+# Mozilla Public License, version 2.0
+
+"""One-step state updates for injection outages using JAX."""
+
+import jax
+import jax.numpy as jnp
+from jax_dataclasses import pytree_dataclass
+from jaxtyping import Array, Complex128, Float, Int
+
+from .lodf_branches import _compute_monitored_branch_currents, _prepare_monitored_branch_pack
+from .lodf_voltages import build_monitor_rows
+
+# ruff: noqa: PLR0913
+
+
+@pytree_dataclass
+class InjectionOutageMonitoredResults:
+    """One-step post-contingency results restricted to monitored elements."""
+
+    n_1_theta: Float[Array, "... n_contingencies n_buses_monitored"]
+    n_1_voltage: Float[Array, "... n_contingencies n_buses_monitored"]
+
+    n_1_p_from: Float[Array, "... n_contingencies n_branches_monitored"]
+    n_1_p_to: Float[Array, "... n_contingencies n_branches_monitored"]
+    n_1_q_from: Float[Array, "... n_contingencies n_branches_monitored"]
+    n_1_q_to: Float[Array, "... n_contingencies n_branches_monitored"]
+    n_1_i_from: Complex128[Array, "... n_contingencies n_branches_monitored"]
+    n_1_i_to: Complex128[Array, "... n_contingencies n_branches_monitored"]
+
+
+def _collect_single_injection_outage_mismatch(
+    outage_injection_indices: Int[jnp.ndarray, " n_outages"],
+    jacobian_size: int,
+    injection_to_bus: Int[jnp.ndarray, " n_injections"],
+    injection_active_power: Float[jnp.ndarray, " n_injections"],
+    injection_reactive_power: Float[jnp.ndarray, " n_injections"],
+    angle_component_indices: Int[jnp.ndarray, " n_buses"],
+    magnitude_component_indices: Int[jnp.ndarray, " n_buses"],
+    dtype: jnp.dtype,
+) -> Float[jnp.ndarray, " n_eq"]:
+    """Assemble the outage-induced mismatch vector for one contingency.
+
+    Parameters
+    ----------
+    outage_injection_indices : Int[jnp.ndarray, " n_outages"]
+            Indices of injections to disconnect for one contingency. Negative entries
+            are ignored and can be used as padding.
+    jacobian_size : int
+            Number of equations in the Jacobian system.
+    injection_to_bus : Int[jnp.ndarray, " n_injections"]
+            Map from injection index to hosting bus index.
+    injection_active_power : Float[jnp.ndarray, " n_injections"]
+            Active power associated with each injection.
+    injection_reactive_power : Float[jnp.ndarray, " n_injections"]
+            Reactive power associated with each injection.
+    angle_component_indices : Int[jnp.ndarray, " n_buses"]
+            Map from bus index to active-power mismatch component index.
+    magnitude_component_indices : Int[jnp.ndarray, " n_buses"]
+            Map from bus index to reactive-power mismatch component index.
+    dtype : jnp.dtype
+            Target dtype of the assembled mismatch vector.
+
+    Returns
+    -------
+    Float[jnp.ndarray, " n_eq"]
+            Mismatch vector in Jacobian component ordering.
+    """
+    valid_outage_mask = outage_injection_indices >= 0
+    safe_outages = jnp.where(valid_outage_mask, outage_injection_indices, 0)
+    buses = injection_to_bus[safe_outages]
+
+    mismatch = jnp.zeros((jacobian_size,), dtype=dtype)
+
+    p_idx = angle_component_indices[buses]
+    valid_p = valid_outage_mask & (p_idx >= 0)
+    safe_p_idx = jnp.where(valid_p, p_idx, 0)
+    p_values = jnp.where(valid_p, -injection_active_power[safe_outages], 0.0)
+    mismatch = mismatch.at[safe_p_idx].add(p_values)
+
+    q_idx = magnitude_component_indices[buses]
+    valid_q = valid_outage_mask & (q_idx >= 0)
+    safe_q_idx = jnp.where(valid_q, q_idx, 0)
+    q_values = jnp.where(valid_q, -injection_reactive_power[safe_outages], 0.0)
+    mismatch = mismatch.at[safe_q_idx].add(q_values)
+
+    return mismatch
+
+
+def _collect_injection_outage_mismatch(
+    jacobian_size: int,
+    outage_injection_indices: Int[jnp.ndarray, " n_contingencies n_outages"],
+    injection_to_bus: Int[jnp.ndarray, " n_injections"],
+    injection_active_power: Float[jnp.ndarray, " n_injections"],
+    injection_reactive_power: Float[jnp.ndarray, " n_injections"],
+    angle_component_indices: Int[jnp.ndarray, " n_buses"],
+    magnitude_component_indices: Int[jnp.ndarray, " n_buses"],
+    dtype: jnp.dtype,
+) -> Float[jnp.ndarray, " n_contingencies n_eq"]:
+    """Assemble outage-induced mismatch vectors in Jacobian ordering.
+
+    Parameters
+    ----------
+    jacobian_size : int
+            Number of equations in the Jacobian system.
+    outage_injection_indices : Int[jnp.ndarray, " n_contingencies n_outages"]
+            Batched injection outage indices. Each row represents one contingency and
+            may include negative padding entries.
+    injection_to_bus : Int[jnp.ndarray, " n_injections"]
+            Map from injection index to hosting bus index.
+    injection_active_power : Float[jnp.ndarray, " n_injections"]
+            Active power associated with each injection.
+    injection_reactive_power : Float[jnp.ndarray, " n_injections"]
+            Reactive power associated with each injection.
+    angle_component_indices : Int[jnp.ndarray, " n_buses"]
+            Map from bus index to active-power mismatch component index.
+    magnitude_component_indices : Int[jnp.ndarray, " n_buses"]
+            Map from bus index to reactive-power mismatch component index.
+    dtype : jnp.dtype
+            Target dtype of the assembled mismatch matrix.
+
+    Returns
+    -------
+    Float[jnp.ndarray, " n_contingencies n_eq"]
+            Batched mismatch vectors in Jacobian component ordering.
+    """
+
+    def collect_single(outage_row: Int[jnp.ndarray, " n_outages"]) -> Float[jnp.ndarray, " n_eq"]:
+        return _collect_single_injection_outage_mismatch(
+            outage_injection_indices=outage_row,
+            jacobian_size=jacobian_size,
+            injection_to_bus=injection_to_bus,
+            injection_active_power=injection_active_power,
+            injection_reactive_power=injection_reactive_power,
+            angle_component_indices=angle_component_indices,
+            magnitude_component_indices=magnitude_component_indices,
+            dtype=dtype,
+        )
+
+    return jax.vmap(collect_single)(outage_injection_indices)
+
+
+@jax.jit
+def _calculate_injection_outage_dx(
+    jacobian_inv_transposed: Float[jnp.ndarray, " n_eq n_eq"],
+    mismatch: Float[jnp.ndarray, " n_contingencies n_eq"],
+) -> Float[jnp.ndarray, " n_contingencies n_eq"]:
+    """Map outage mismatch vectors to one-step state increments.
+
+    Parameters
+    ----------
+    jacobian_inv_transposed : Float[jnp.ndarray, " n_eq n_eq"]
+            Transposed inverse Jacobian at the base operating point.
+    mismatch : Float[jnp.ndarray, " n_contingencies n_eq"]
+            Batched outage-induced mismatch vectors.
+
+    Returns
+    -------
+    Float[jnp.ndarray, " n_contingencies n_eq"]
+            Batched one-step state increments in Jacobian ordering.
+    """
+    return -(mismatch @ jacobian_inv_transposed)
+
+
+@jax.jit
+def _calculate_monitor_bus_state_updates(
+    jacobian_inv_transposed: Float[jnp.ndarray, " n_eq n_eq"],
+    mismatch: Float[jnp.ndarray, " n_contingencies n_eq"],
+    theta_rows: Int[jnp.ndarray, " n_mon_bus"],
+    vm_rows: Int[jnp.ndarray, " n_mon_bus"],
+    theta_mask: Float[jnp.ndarray, " n_mon_bus"],
+    vm_mask: Float[jnp.ndarray, " n_mon_bus"],
+) -> tuple[
+    Float[jnp.ndarray, " n_contingencies n_mon_bus"],
+    Float[jnp.ndarray, " n_contingencies n_mon_bus"],
+]:
+    """Map outage mismatch vectors to monitored bus state increments only.
+
+    Parameters
+    ----------
+    jacobian_inv_transposed : Float[jnp.ndarray, " n_eq n_eq"]
+            Transposed inverse Jacobian at the base operating point.
+    mismatch : Float[jnp.ndarray, " n_contingencies n_eq"]
+            Batched outage-induced mismatch vectors.
+    theta_rows : Int[jnp.ndarray, " n_mon_bus"]
+            Safe Jacobian row indices for monitored bus angles.
+    vm_rows : Int[jnp.ndarray, " n_mon_bus"]
+            Safe Jacobian row indices for monitored bus magnitudes.
+    theta_mask : Float[jnp.ndarray, " n_mon_bus"]
+            Mask for monitored buses with angle states in the Jacobian.
+    vm_mask : Float[jnp.ndarray, " n_mon_bus"]
+            Mask for monitored buses with magnitude states in the Jacobian.
+
+    Returns
+    -------
+    tuple[Float[jnp.ndarray, " n_contingencies n_mon_bus"], Float[jnp.ndarray, " n_contingencies n_mon_bus"]]
+            Monitored angle and magnitude increments for each contingency.
+    """
+    dtype = jacobian_inv_transposed.dtype
+    theta_mask_d = theta_mask.astype(dtype)
+    vm_mask_d = vm_mask.astype(dtype)
+    theta_dx = -(mismatch @ jacobian_inv_transposed[:, theta_rows]) * theta_mask_d[None, :]
+    vm_dx = -(mismatch @ jacobian_inv_transposed[:, vm_rows]) * vm_mask_d[None, :]
+    return theta_dx, vm_dx
+
+
+def non_voltage_regulating_injection_outage_monitor_buses(
+    jacobian_inv_transposed: Float[jnp.ndarray, " n_eq n_eq"],
+    outage_injection_indices: Int[jnp.ndarray, " n_contingencies n_outages"],
+    injection_to_bus: Int[jnp.ndarray, " n_injections"],
+    injection_active_power: Float[jnp.ndarray, " n_injections"],
+    injection_reactive_power: Float[jnp.ndarray, " n_injections"],
+    angle_component_indices: Int[jnp.ndarray, " n_buses"],
+    magnitude_component_indices: Int[jnp.ndarray, " n_buses"],
+    monitor_bus_indices: Int[jnp.ndarray, " n_mon_bus"],
+    v_mag_hat: Float[jnp.ndarray, " n_buses"],
+    theta_hat: Float[jnp.ndarray, " n_buses"],
+) -> tuple[
+    Float[jnp.ndarray, " n_contingencies n_mon_bus"],
+    Float[jnp.ndarray, " n_contingencies n_mon_bus"],
+]:
+    """Compute monitored post-contingency bus voltages for injection outages only.
+
+    Parameters
+    ----------
+    jacobian_inv_transposed : Float[jnp.ndarray, " n_eq n_eq"]
+            Base transposed inverse Jacobian at the hot-start operating point.
+    outage_injection_indices : Int[jnp.ndarray, " n_contingencies n_outages"]
+            Batched indices of injections to disconnect. Each row is one contingency.
+    injection_to_bus : Int[jnp.ndarray, " n_injections"]
+            Map from injection index to hosting bus index.
+    injection_active_power : Float[jnp.ndarray, " n_injections"]
+            Active power of each injection at the hot-start operating point.
+    injection_reactive_power : Float[jnp.ndarray, " n_injections"]
+            Reactive power of each injection at the hot-start operating point.
+    angle_component_indices : Int[jnp.ndarray, " n_buses"]
+            Map from bus index to angle component index in the Jacobian. -1 if not present.
+    magnitude_component_indices : Int[jnp.ndarray, " n_buses"]
+            Map from bus index to magnitude component index in the Jacobian. -1 if not present.
+    monitor_bus_indices : Int[jnp.ndarray, " n_mon_bus"]
+            Indices of buses whose post-contingency voltages should be returned.
+    v_mag_hat : Float[jnp.ndarray, " n_buses"]
+            Base-case voltage magnitudes.
+    theta_hat : Float[jnp.ndarray, " n_buses"]
+            Base-case voltage angles.
+
+    Returns
+    -------
+    tuple[Float[jnp.ndarray, " n_contingencies n_mon_bus"], Float[jnp.ndarray, " n_contingencies n_mon_bus"]]
+            Post-contingency monitored voltage angles and magnitudes.
+    """
+    dtype = jacobian_inv_transposed.dtype
+    mismatch = _collect_injection_outage_mismatch(
+        jacobian_size=jacobian_inv_transposed.shape[0],
+        outage_injection_indices=outage_injection_indices,
+        injection_to_bus=injection_to_bus,
+        injection_active_power=injection_active_power,
+        injection_reactive_power=injection_reactive_power,
+        angle_component_indices=angle_component_indices,
+        magnitude_component_indices=magnitude_component_indices,
+        dtype=dtype,
+    )
+    theta_rows, vm_rows, theta_mask, vm_mask = build_monitor_rows(
+        angle_component_indices=angle_component_indices,
+        magnitude_component_indices=magnitude_component_indices,
+        monitor_bus_indices=monitor_bus_indices,
+    )
+    theta_dx, vm_dx = _calculate_monitor_bus_state_updates(
+        jacobian_inv_transposed=jacobian_inv_transposed,
+        mismatch=mismatch,
+        theta_rows=theta_rows,
+        vm_rows=vm_rows,
+        theta_mask=theta_mask,
+        vm_mask=vm_mask,
+    )
+    base_theta0 = jnp.take(theta_hat, monitor_bus_indices, axis=0).astype(dtype)
+    base_vm0 = jnp.take(v_mag_hat, monitor_bus_indices, axis=0).astype(dtype)
+    return base_theta0[None, :] + theta_dx, base_vm0[None, :] + vm_dx
+
+
+def non_voltage_regulating_injection_outage_monitored(
+    jacobian_inv_transposed: Float[jnp.ndarray, " n_eq n_eq"],
+    outage_injection_indices: Int[jnp.ndarray, " n_contingencies n_outages"],
+    injection_to_bus: Int[jnp.ndarray, " n_injections"],
+    injection_active_power: Float[jnp.ndarray, " n_injections"],
+    injection_reactive_power: Float[jnp.ndarray, " n_injections"],
+    angle_component_indices: Int[jnp.ndarray, " n_buses"],
+    magnitude_component_indices: Int[jnp.ndarray, " n_buses"],
+    monitor_bus_indices: Int[jnp.ndarray, " n_mon_bus"],
+    v_mag_hat: Float[jnp.ndarray, " n_buses"],
+    theta_hat: Float[jnp.ndarray, " n_buses"],
+    branch_from: Int[jnp.ndarray, " n_branches"],
+    branch_to: Int[jnp.ndarray, " n_branches"],
+    y_ff: Complex128[jnp.ndarray, " n_branches"],
+    y_ft: Complex128[jnp.ndarray, " n_branches"],
+    y_tf: Complex128[jnp.ndarray, " n_branches"],
+    y_tt: Complex128[jnp.ndarray, " n_branches"],
+    monitor_branch_indices: Int[jnp.ndarray, " n_mon_br"],
+    bus_to_mon_index: Int[jnp.ndarray, " n_buses"],
+) -> InjectionOutageMonitoredResults:
+    """Compute monitored bus states and branch flows for injection outages.
+
+    Parameters
+    ----------
+    jacobian_inv_transposed : Float[jnp.ndarray, " n_eq n_eq"]
+            Base transposed inverse Jacobian at the hot-start operating point.
+    outage_injection_indices : Int[jnp.ndarray, " n_contingencies n_outages"]
+            Batched indices of injections to disconnect. Each row is one contingency.
+    injection_to_bus : Int[jnp.ndarray, " n_injections"]
+            Map from injection index to hosting bus index.
+    injection_active_power : Float[jnp.ndarray, " n_injections"]
+            Active power of each injection at the hot-start operating point.
+    injection_reactive_power : Float[jnp.ndarray, " n_injections"]
+            Reactive power of each injection at the hot-start operating point.
+    angle_component_indices : Int[jnp.ndarray, " n_buses"]
+            Map from bus index to angle component index in the Jacobian. -1 if not present.
+    magnitude_component_indices : Int[jnp.ndarray, " n_buses"]
+            Map from bus index to magnitude component index in the Jacobian. -1 if not present.
+    monitor_bus_indices : Int[jnp.ndarray, " n_mon_bus"]
+            Indices of buses whose post-contingency voltages should be returned.
+    v_mag_hat : Float[jnp.ndarray, " n_buses"]
+            Base-case voltage magnitudes.
+    theta_hat : Float[jnp.ndarray, " n_buses"]
+            Base-case voltage angles.
+    branch_from : Int[jnp.ndarray, " n_branches"]
+            From bus index of each branch.
+    branch_to : Int[jnp.ndarray, " n_branches"]
+            To bus index of each branch.
+    y_ff : Complex128[jnp.ndarray, " n_branches"]
+            From-from branch admittances.
+    y_ft : Complex128[jnp.ndarray, " n_branches"]
+            From-to branch admittances.
+    y_tf : Complex128[jnp.ndarray, " n_branches"]
+            To-from branch admittances.
+    y_tt : Complex128[jnp.ndarray, " n_branches"]
+            To-to branch admittances.
+    monitor_branch_indices : Int[jnp.ndarray, " n_mon_br"]
+            Indices of monitored branches.
+    bus_to_mon_index : Int[jnp.ndarray, " n_buses"]
+            Mapping from global bus index to monitored-bus position. Use -1 for unmonitored buses.
+
+    Returns
+    -------
+    InjectionOutageMonitoredResults
+            Post-contingency monitored bus states and monitored branch flows.
+    """
+    dtype = jacobian_inv_transposed.dtype
+    theta_all, vm_all = non_voltage_regulating_injection_outage_monitor_buses(
+        jacobian_inv_transposed=jacobian_inv_transposed,
+        outage_injection_indices=outage_injection_indices,
+        injection_to_bus=injection_to_bus,
+        injection_active_power=injection_active_power,
+        injection_reactive_power=injection_reactive_power,
+        angle_component_indices=angle_component_indices,
+        magnitude_component_indices=magnitude_component_indices,
+        monitor_bus_indices=monitor_bus_indices,
+        v_mag_hat=v_mag_hat,
+        theta_hat=theta_hat,
+    )
+    (
+        _,
+        y_ff_mon,
+        y_ft_mon,
+        y_tf_mon,
+        y_tt_mon,
+        f_pos_safe,
+        t_pos_safe,
+        end_mask,
+    ) = _prepare_monitored_branch_pack(
+        branch_from=branch_from,
+        branch_to=branch_to,
+        y_ff=y_ff,
+        y_ft=y_ft,
+        y_tf=y_tf,
+        y_tt=y_tt,
+        monitor_branch_indices=monitor_branch_indices,
+        bus_to_mon_index=bus_to_mon_index,
+        dtype=dtype,
+    )
+    v_from_all, v_to_all, i_from_all, i_to_all = _compute_monitored_branch_currents(
+        theta_all=theta_all,
+        vm_all=vm_all,
+        f_pos_safe=f_pos_safe,
+        t_pos_safe=t_pos_safe,
+        y_ff_mon=y_ff_mon,
+        y_ft_mon=y_ft_mon,
+        y_tf_mon=y_tf_mon,
+        y_tt_mon=y_tt_mon,
+        dtype=dtype,
+    )
+    complex_dtype = y_ff_mon.dtype
+    end_mask_complex = end_mask.astype(complex_dtype)[None, :]
+    end_mask_real = end_mask.astype(dtype)[None, :]
+    s_from_all = v_from_all * jnp.conj(i_from_all) * end_mask_complex
+    s_to_all = v_to_all * jnp.conj(i_to_all) * end_mask_complex
+    return InjectionOutageMonitoredResults(
+        n_1_theta=theta_all,
+        n_1_voltage=vm_all,
+        n_1_p_from=s_from_all.real.astype(dtype) * end_mask_real,
+        n_1_p_to=s_to_all.real.astype(dtype) * end_mask_real,
+        n_1_q_from=s_from_all.imag.astype(dtype) * end_mask_real,
+        n_1_q_to=s_to_all.imag.astype(dtype) * end_mask_real,
+        n_1_i_from=i_from_all * end_mask_complex,
+        n_1_i_to=i_to_all * end_mask_complex,
+    )
+
+
+def non_voltage_regulating_injection_outage_dx(
+    jacobian_inv_transposed: Float[jnp.ndarray, " n_eq n_eq"],
+    outage_injection_indices: Int[jnp.ndarray, " n_contingencies n_outages"],
+    injection_to_bus: Int[jnp.ndarray, " n_injections"],
+    injection_active_power: Float[jnp.ndarray, " n_injections"],
+    injection_reactive_power: Float[jnp.ndarray, " n_injections"],
+    angle_component_indices: Int[jnp.ndarray, " n_buses"],
+    magnitude_component_indices: Int[jnp.ndarray, " n_buses"],
+) -> Float[jnp.ndarray, " n_contingencies n_eq"]:
+    """One-step state increment caused by disconnecting non voltage regulating injections.
+
+    The Jacobian is unchanged for a pure injection outage because network
+    topology and admittances remain fixed. The outage only changes the nodal
+    power mismatch at the buses hosting the disconnected injections.
+
+    This helper assumes the bus type pattern used by the Jacobian remains
+    unchanged across the outage. That is appropriate for load outages and other
+    contingencies that do not alter voltage-control structure.
+
+    Note: this function still works for voltage regulating injection outages,
+    but the PV->PQ bus type change is not accounted for.
+
+    Parameters
+    ----------
+    jacobian_inv_transposed : Float[jnp.ndarray, " n_eq n_eq"]
+            Base transposed inverse Jacobian at the hot-start operating point.
+    outage_injection_indices : Int[jnp.ndarray, " n_contingencies n_outages"]
+            Batched indices of injections to disconnect. Each row is one contingency.
+    injection_to_bus : Int[jnp.ndarray, " n_injections"]
+            Map from injection index to hosting bus index.
+    injection_active_power : Float[jnp.ndarray, " n_injections"]
+            Active power of each injection at the hot-start operating point.
+    injection_reactive_power : Float[jnp.ndarray, " n_injections"]
+            Reactive power of each injection at the hot-start operating point.
+    angle_component_indices : Int[jnp.ndarray, " n_buses"]
+            Map from bus index to angle component index in the Jacobian. -1 if not present.
+    magnitude_component_indices : Int[jnp.ndarray, " n_buses"]
+            Map from bus index to magnitude component index in the Jacobian. -1 if not present.
+
+    Returns
+    -------
+    Float[jnp.ndarray, " n_contingencies n_eq"]
+            One-step state increment for each outage in Jacobian ordering.
+    """
+    dtype = jacobian_inv_transposed.dtype
+    mismatch = _collect_injection_outage_mismatch(
+        jacobian_size=jacobian_inv_transposed.shape[0],
+        outage_injection_indices=outage_injection_indices,
+        injection_to_bus=injection_to_bus,
+        injection_active_power=injection_active_power,
+        injection_reactive_power=injection_reactive_power,
+        angle_component_indices=angle_component_indices,
+        magnitude_component_indices=magnitude_component_indices,
+        dtype=dtype,
+    )
+    return _calculate_injection_outage_dx(jacobian_inv_transposed=jacobian_inv_transposed, mismatch=mismatch)
diff --git a/src/dc_plus/jax/low_rank_helper.py b/src/dc_plus/jax/low_rank_helper.py
index 7089ffe..85f1af5 100644
--- a/src/dc_plus/jax/low_rank_helper.py
+++ b/src/dc_plus/jax/low_rank_helper.py
@@ -128,7 +128,9 @@ def _compute_branch_delta_submatrix_from_admittance(
     dqt_dvt = -2.0 * v_t * b_tt + v_f * (g_tf * sin_tf - b_tf * cos_tf)
     dqt_dvf = v_t * (g_tf * sin_tf - b_tf * cos_tf)
 
-    dtype = jnp.result_type(v_mag_from, v_mag_to, theta_from, theta_to, y_ff)
+    # Jacobian block entries are real-valued; avoid propagating complex dtypes
+    # from admittance inputs into downstream scatter-add updates.
+    dtype_input = jnp.result_type(v_mag_from, v_mag_to, theta_from, theta_to, g_ff)
 
     delta = jnp.array(
         [
@@ -137,7 +139,7 @@ def _compute_branch_delta_submatrix_from_admittance(
             [-dqf_dthf, -dqf_dtht, -dqf_dvf, -dqf_dvt],
             [-dqt_dthf, -dqt_dtht, -dqt_dvf, -dqt_dvt],
         ],
-        dtype=dtype,
+        dtype=dtype_input,
     )
     return delta
 
diff --git a/src/dc_plus/numpy/__init__.py b/src/dc_plus/numpy/__init__.py
index 9c761f9..dc4b87d 100644
--- a/src/dc_plus/numpy/__init__.py
+++ b/src/dc_plus/numpy/__init__.py
@@ -10,3 +10,7 @@
 This implementation is intended for reference and testing purposes only.
 It is currently not optimized for performance, but for readability.
 """
+
+from .injection_outage import non_voltage_regulating_injection_outage_dx
+
+__all__ = ["non_voltage_regulating_injection_outage_dx"]
diff --git a/src/dc_plus/numpy/injection_outage.py b/src/dc_plus/numpy/injection_outage.py
new file mode 100644
index 0000000..666eaa2
--- /dev/null
+++ b/src/dc_plus/numpy/injection_outage.py
@@ -0,0 +1,170 @@
+# Copyright 2026 50Hertz Transmission GmbH and Elia Transmission Belgium SA/NV
+#
+# This Source Code Form is subject to the terms of the Mozilla Public License, v. 2.0.
+# If a copy of the MPL was not distributed with this file,
+# you can obtain one at https://mozilla.org/MPL/2.0/.
+# Mozilla Public License, version 2.0
+
+"""One-step state updates for injection outages."""
+
+import numpy as np
+from jaxtyping import Float, Int
+
+
+def _collect_injection_outage_mismatch(
+    jacobian_size: int,
+    outage_injection_indices: Int[np.ndarray, " n_contingencies n_outages"],
+    injection_to_bus: Int[np.ndarray, " n_injections"],
+    injection_active_power: Float[np.ndarray, " n_injections"],
+    injection_reactive_power: Float[np.ndarray, " n_injections"],
+    angle_component_indices: Int[np.ndarray, " n_buses"],
+    magnitude_component_indices: Int[np.ndarray, " n_buses"],
+    dtype: np.dtype,
+) -> Float[np.ndarray, " n_contingencies n_eq"]:
+    """Assemble outage-induced mismatch vectors in Jacobian ordering.
+
+    Parameters
+    ----------
+    jacobian_size : int
+        Number of equations in the Jacobian system.
+    outage_injection_indices : Int[np.ndarray, " n_contingencies n_outages"]
+        Batched injection outage indices. Each row represents one contingency and
+        may include negative padding entries.
+    injection_to_bus : Int[np.ndarray, " n_injections"]
+        Map from injection index to hosting bus index.
+    injection_active_power : Float[np.ndarray, " n_injections"]
+        Active power associated with each injection.
+    injection_reactive_power : Float[np.ndarray, " n_injections"]
+        Reactive power associated with each injection.
+    angle_component_indices : Int[np.ndarray, " n_buses"]
+        Map from bus index to active-power mismatch component index.
+    magnitude_component_indices : Int[np.ndarray, " n_buses"]
+        Map from bus index to reactive-power mismatch component index.
+    dtype : np.dtype
+        Target dtype of the assembled mismatch matrix.
+
+    Returns
+    -------
+    Float[np.ndarray, " n_contingencies n_eq"]
+        Batched mismatch vectors in Jacobian component ordering.
+    """
+    outage_indices = outage_injection_indices
+    mismatch = np.zeros((outage_indices.shape[0], jacobian_size), dtype=dtype)
+    if outage_indices.size == 0:
+        return mismatch
+
+    valid_outage_mask = outage_indices >= 0
+    if not np.any(valid_outage_mask):
+        return mismatch
+
+    valid_rows, valid_cols = np.nonzero(valid_outage_mask)
+    valid_outages = outage_indices[valid_rows, valid_cols]
+    if np.any(valid_outages >= injection_to_bus.size):
+        raise IndexError("Injection outage index is out of bounds")
+
+    buses = injection_to_bus[valid_outages]
+
+    p_idx = angle_component_indices[buses]
+    valid_p = p_idx >= 0
+    if np.any(valid_p):
+        np.add.at(
+            mismatch,
+            (valid_rows[valid_p], p_idx[valid_p]),
+            -injection_active_power[valid_outages][valid_p],
+        )
+
+    q_idx = magnitude_component_indices[buses]
+    valid_q = q_idx >= 0
+    if np.any(valid_q):
+        np.add.at(
+            mismatch,
+            (valid_rows[valid_q], q_idx[valid_q]),
+            -injection_reactive_power[valid_outages][valid_q],
+        )
+
+    return mismatch
+
+
+def _calculate_injection_outage_dx(
+    jacobian_inv: Float[np.ndarray, " n_eq n_eq"],
+    mismatch: Float[np.ndarray, " n_contingencies n_eq"],
+) -> Float[np.ndarray, " n_contingencies n_eq"]:
+    """Map outage mismatch vectors to one-step state increments.
+
+    Parameters
+    ----------
+    jacobian_inv : Float[np.ndarray, " n_eq n_eq"]
+        Inverse Jacobian at the base operating point.
+    mismatch : Float[np.ndarray, " n_contingencies n_eq"]
+        Batched outage-induced mismatch vectors.
+
+    Returns
+    -------
+    Float[np.ndarray, " n_contingencies n_eq"]
+        Batched one-step state increments in Jacobian ordering.
+    """
+    return -(mismatch @ jacobian_inv.T)
+
+
+def non_voltage_regulating_injection_outage_dx(
+    jacobian_inv: Float[np.ndarray, " n_eq n_eq"],
+    outage_injection_indices: Int[np.ndarray, " n_contingencies n_outages"],
+    injection_to_bus: Int[np.ndarray, " n_injections"],
+    injection_active_power: Float[np.ndarray, " n_injections"],
+    injection_reactive_power: Float[np.ndarray, " n_injections"],
+    angle_component_indices: Int[np.ndarray, " n_buses"],
+    magnitude_component_indices: Int[np.ndarray, " n_buses"],
+) -> Float[np.ndarray, " n_contingencies n_eq"]:
+    """One-step state increment caused by disconnecting non voltage regulating injections.
+
+    The Jacobian is unchanged for a pure injection outage because network
+    topology and admittances remain fixed. The outage only changes the nodal
+    power mismatch at the buses hosting the disconnected injections.
+
+    This helper assumes the bus type pattern used by the Jacobian remains
+    unchanged across the outage. That is appropriate for load outages and other
+    contingencies that do not alter voltage-control structure.
+
+    Note: this function still works for voltage regulating injection outages,
+    but the PV->PQ bus type change is not accounted for.
+
+    Parameters
+    ----------
+    jacobian_inv : Float[np.ndarray, " n_eq n_eq"]
+        Base inverse Jacobian at the hot-start operating point.
+    outage_injection_indices : Int[np.ndarray, " n_contingencies n_outages"]
+        Batched indices of injections to disconnect. Each row is one
+        contingency and may contain multiple outages. Negative entries are
+        ignored and can be used as padding for rows with different outage
+        counts.
+    injection_to_bus : Int[np.ndarray, " n_injections"]
+        Mapping from injection index to bus index.
+    injection_active_power : Float[np.ndarray, " n_injections"]
+        Active power removed when each injection is disconnected.
+        In practice this is typically the injection setpoint, or the solved
+        injection power when no separate setpoint is available.
+    injection_reactive_power : Float[np.ndarray, " n_injections"]
+        Reactive power removed when each injection is disconnected.
+        In practice this is typically the injection setpoint, or the solved
+        injection power when no separate setpoint is available.
+    angle_component_indices : Int[np.ndarray, " n_buses"]
+        Bus to active-power Jacobian row mapping.
+    magnitude_component_indices : Int[np.ndarray, " n_buses"]
+        Bus to reactive-power Jacobian row mapping.
+
+    Returns
+    -------
+    Float[np.ndarray, " n_contingencies n_eq"]
+        Batched one-step update vectors in Jacobian ordering.
+    """
+    mismatch = _collect_injection_outage_mismatch(
+        jacobian_size=jacobian_inv.shape[0],
+        outage_injection_indices=outage_injection_indices,
+        injection_to_bus=injection_to_bus,
+        injection_active_power=injection_active_power,
+        injection_reactive_power=injection_reactive_power,
+        angle_component_indices=angle_component_indices,
+        magnitude_component_indices=magnitude_component_indices,
+        dtype=jacobian_inv.dtype,
+    )
+    return _calculate_injection_outage_dx(jacobian_inv=jacobian_inv, mismatch=mismatch)
diff --git a/src/dc_plus/numpy/low_rank_helper.py b/src/dc_plus/numpy/low_rank_helper.py
index 9299b4d..03addd7 100644
--- a/src/dc_plus/numpy/low_rank_helper.py
+++ b/src/dc_plus/numpy/low_rank_helper.py
@@ -209,7 +209,7 @@ def _compute_branch_delta_submatrix_from_admittance(
     dqt_dvt = -2.0 * v_t * b_tt + v_f * (g_tf * sin_tf - b_tf * cos_tf)
     dqt_dvf = v_t * (g_tf * sin_tf - b_tf * cos_tf)
 
-    dtype = np.result_type(v_mag_from, v_mag_to, theta_from, theta_to, y_ff)
+    dtype = np.result_type(v_mag_from, v_mag_to, theta_from, theta_to, g_ff)
 
     delta = np.array(
         [
diff --git a/src/dc_plus/preprocess/create_network_data.py b/src/dc_plus/preprocess/create_network_data.py
index d7f4470..9efafd0 100644
--- a/src/dc_plus/preprocess/create_network_data.py
+++ b/src/dc_plus/preprocess/create_network_data.py
@@ -140,6 +140,7 @@ def _create_network_data(
         bus_active_power=bus_active_power,
         bus_reactive_power=bus_reactive_power,
         bus_type=buses["bus_type"].values.astype(int),
+        bus_is_angle_reference=buses["is_angle_reference"].values.astype(bool),
         injection_to_bus=injections["bus_index"].values.astype(int),
         injection_active_power=injections["p"].values,
         injection_reactive_power=injections["q"].values,
@@ -187,7 +188,12 @@ def create_network_data(
     injections = _get_injections_powsybl(network)
     shunts = _get_shunts_powsybl(network)
     slack_id = network.get_extensions("slackTerminal")["bus_id"].values[0]
-    buses = _get_buses_powsybl(net=network, slack_id=slack_id, injections=injections)
+    references = network.get_extensions("referencePriorities")
+    if len(references) > 0:
+        reference_id = network.get_extensions("referencePriorities").index[0]
+    else:
+        reference_id = slack_id
+    buses = _get_buses_powsybl(net=network, slack_id=slack_id, injections=injections, reference_id=reference_id)
     limits = _get_limits_parameter_powsybl(network)
 
     return _create_network_data(
diff --git a/src/dc_plus/preprocess/preprocess_jacobian_bsdf.py b/src/dc_plus/preprocess/preprocess_jacobian_bsdf.py
index dbf7411..5f7c5de 100644
--- a/src/dc_plus/preprocess/preprocess_jacobian_bsdf.py
+++ b/src/dc_plus/preprocess/preprocess_jacobian_bsdf.py
@@ -302,6 +302,12 @@ def _pad_axis0(array: np.ndarray, value: float) -> np.ndarray:
             np.full(split_count, BusType.PQ, dtype=dynamic_network_data.bus_type.dtype),
         )
     )
+    extended_bus_is_angle_reference = np.concatenate(
+        (
+            dynamic_network_data.bus_is_angle_reference,
+            np.zeros(split_count, dtype=dynamic_network_data.bus_is_angle_reference.dtype),
+        )
+    )
 
     return replace(
         dynamic_network_data,
@@ -310,6 +316,7 @@ def _pad_axis0(array: np.ndarray, value: float) -> np.ndarray:
         bus_active_power=extended_bus_active_power,
         bus_reactive_power=extended_bus_reactive_power,
         bus_type=extended_bus_type,
+        bus_is_angle_reference=extended_bus_is_angle_reference,
     )
 
 
diff --git a/tests/importing/powsybl/test_powsybl_import.py b/tests/importing/powsybl/test_powsybl_import.py
index 3d7df2c..f8cfa51 100644
--- a/tests/importing/powsybl/test_powsybl_import.py
+++ b/tests/importing/powsybl/test_powsybl_import.py
@@ -39,6 +39,8 @@
     _get_trafo_parameter,
 )
 
+from dc_plus.interfaces.network_information import BusType
+
 
 def test_get_tie_line_parameter():
     net = pypowsybl.network.create_eurostag_tutorial_example1_with_tie_lines_and_areas()
@@ -188,17 +190,18 @@ def test_get_buses():
     net = create_complex_grid_battery_hvdc_svc_3w_trafo()
     injections = _get_injections_powsybl(net)
     slack_id = net.get_extensions("slackTerminal")["bus_id"].values[0]
-    buses = _get_buses_powsybl(net=net, slack_id=slack_id, injections=injections)
+    buses = _get_buses_powsybl(net=net, slack_id=slack_id, injections=injections, reference_id=slack_id)
     BusParamSchema.validate(buses)
     bus_powsybl = net.get_buses(attributes=[])
     dangling_count = len(_get_dangling_bus_ids(net))
     assert len(bus_powsybl) + dangling_count == len(buses)
-    assert 0 in list(buses["bus_type"].values)  # slack bus
-    assert 1 in list(buses["bus_type"].values)  # PV bus
-    assert 2 in list(buses["bus_type"].values)  # PQ bus
-    assert np.sum(buses["bus_type"] == 0) == 1  # only one slack bus
-    assert np.sum(buses["bus_type"] == 1) == 2  # there are two PV buses
-    assert np.sum(buses["bus_type"] == 2) >= 1  # at least one PQ bus
+    assert BusType.SLACK in list(buses["bus_type"].values)  # slack bus
+    assert BusType.PV in list(buses["bus_type"].values)  # PV bus
+    assert BusType.PQ in list(buses["bus_type"].values)  # PQ bus
+    assert np.sum(buses["bus_type"] == BusType.SLACK) == 1  # only one slack bus
+    assert np.sum(buses["is_angle_reference"]) == 1  # only one reference bus
+    assert np.sum(buses["bus_type"] == BusType.PV) == 2  # there are two PV buses
+    assert np.sum(buses["bus_type"] == BusType.PQ) >= 1  # at least one PQ bus
 
 
 def test_no_synchronous_components_returns_connected_ids():
diff --git a/tests/importing/test_import_helpers.py b/tests/importing/test_import_helpers.py
index 1e43e76..5633238 100644
--- a/tests/importing/test_import_helpers.py
+++ b/tests/importing/test_import_helpers.py
@@ -23,6 +23,7 @@ def test_remove_isolated_buses_injections():
             "voltage_magnitude": [1.0, 0.98, 1.02],
             "voltage_angle": [0.0, -0.05, 0.03],
             "bus_type": [0, 1, 2],
+            "is_angle_reference": [True, False, False],
             "grid_island_id": [0, 0, 1],
         }
     )
@@ -73,6 +74,7 @@ def test_remove_isolated_branches():
             "voltage_magnitude": [1.0, 1.0, 1.0],
             "voltage_angle": [0.0, 0.0, 0.0],
             "bus_type": [0, 1, 2],
+            "is_angle_reference": [True, False, False],
             "grid_island_id": [0, 0, 1],
         }
     )
@@ -125,6 +127,7 @@ def test_remove_isolated_buses():
             "voltage_magnitude": [1.0, 0.97, 1.01, 0.95],
             "voltage_angle": [0.0, -0.03, 0.02, -0.05],
             "bus_type": [0, 1, 2, 2],
+            "is_angle_reference": [True, False, False, False],
             "grid_island_id": [0, 1, 0, 2],
         }
     )
diff --git a/tests/importing/test_powsybl_vs_pandapower_dynamic.py b/tests/importing/test_powsybl_vs_pandapower_dynamic.py
index e331c58..3bcd137 100644
--- a/tests/importing/test_powsybl_vs_pandapower_dynamic.py
+++ b/tests/importing/test_powsybl_vs_pandapower_dynamic.py
@@ -296,6 +296,7 @@ def _prepare_powsybl_network(
         bus_active_power=bus_active_power,
         bus_reactive_power=bus_reactive_power,
         bus_type=bus_type,
+        bus_is_angle_reference=dynamic_psb.bus_is_angle_reference[bus_perm].copy(),
         injection_to_bus=injection_to_bus,
         injection_active_power=injection_active_power,
         injection_reactive_power=injection_reactive_power,
@@ -384,6 +385,7 @@ def _prepare_pandapower_network(
         bus_active_power=dynamic_pp.bus_active_power[bus_perm].copy(),
         bus_reactive_power=dynamic_pp.bus_reactive_power[bus_perm].copy(),
         bus_type=dynamic_pp.bus_type[bus_perm].copy(),
+        bus_is_angle_reference=dynamic_pp.bus_is_angle_reference[bus_perm].copy(),
         injection_to_bus=dynamic_pp.injection_to_bus[injection_perm].copy(),
         injection_active_power=dynamic_pp.injection_active_power[injection_perm].copy(),
         injection_reactive_power=dynamic_pp.injection_reactive_power[injection_perm].copy(),
diff --git a/tests/jax/test_injection_outage_jax.py b/tests/jax/test_injection_outage_jax.py
new file mode 100644
index 0000000..cc0a021
--- /dev/null
+++ b/tests/jax/test_injection_outage_jax.py
@@ -0,0 +1,241 @@
+# Copyright 2026 50Hertz Transmission GmbH and Elia Transmission Belgium SA/NV
+#
+# This Source Code Form is subject to the terms of the Mozilla Public License, v. 2.0.
+# If a copy of the MPL was not distributed with this file,
+# you can obtain one at https://mozilla.org/MPL/2.0/.
+# Mozilla Public License, version 2.0
+
+import jax
+import jax.numpy as jnp
+import numpy as np
+import pytest
+
+jax.config.update("jax_enable_x64", True)
+jax.config.update("jax_platform_name", "cpu")
+
+from dc_plus.example_grids.pypowsbl.example_grids import PANDAPOWER_NETWORKS_FOR_POWSYBL, POWSYBL_NETWORKS
+from dc_plus.importing.powsybl.powsybl_import import _get_injections_powsybl
+from dc_plus.importing.powsybl.powsybl_network_helpers import _load_test_grid
+from dc_plus.jax.injection_outage import (
+    non_voltage_regulating_injection_outage_monitor_buses,
+    non_voltage_regulating_injection_outage_monitored,
+)
+from dc_plus.jax.injection_outage import non_voltage_regulating_injection_outage_dx as jax_injection_outage_dx
+from dc_plus.numpy.injection_outage import non_voltage_regulating_injection_outage_dx as numpy_injection_outage_dx
+
+powsybl_networks = POWSYBL_NETWORKS
+pandapower_networks = PANDAPOWER_NETWORKS_FOR_POWSYBL
+TOL = 1e-10
+
+
+def _connected_injection_indices(injections):
+    candidates = np.flatnonzero(injections["connected"].to_numpy(dtype=bool))
+    if candidates.size == 0:
+        pytest.skip("No connected injections available for an injection outage test")
+    return candidates
+
+
+def _supported_injection_outage_rows(net, injections, max_rows: int | None = None) -> list[np.ndarray]:
+    injection_powsybl_type = net.get_injections(attributes=["type"])
+    outage_rows: list[np.ndarray] = []
+    for injection_idx in _connected_injection_indices(injections):
+        injection_id = injections.loc[injection_idx, "id_str"]
+        if (injection_powsybl_type.loc[injection_id, "type"]) in (["HVDC_CONVERTER_STATION", "DANGLING_LINE"]):
+            continue
+        outage_rows.append(np.array([injection_idx], dtype=np.int64))
+        if max_rows is not None and len(outage_rows) >= max_rows:
+            break
+    return outage_rows
+
+
+def _pad_outage_rows(outage_rows: list[np.ndarray]) -> np.ndarray:
+    if not outage_rows:
+        return np.empty((0, 0), dtype=np.int64)
+
+    max_outages = max(outage_row.size for outage_row in outage_rows)
+    padded = np.full((len(outage_rows), max_outages), -1, dtype=np.int64)
+    for row_idx, outage_row in enumerate(outage_rows):
+        padded[row_idx, : outage_row.size] = outage_row
+    return padded
+
+
+def test_non_voltage_regulating_injection_outage_jax_batch_multi_outage() -> None:
+    jacobian_inv = np.eye(4, dtype=float)
+    outage_batch = np.array([[0, 2], [1, -1], [-1, -1]], dtype=np.int64)
+    injection_to_bus = np.array([0, 1, 0], dtype=np.int64)
+    injection_active_power = np.array([1.0, 2.0, 3.0], dtype=float)
+    injection_reactive_power = np.array([10.0, 20.0, 30.0], dtype=float)
+    angle_component_indices = np.array([0, 1], dtype=np.int64)
+    magnitude_component_indices = np.array([2, 3], dtype=np.int64)
+
+    dx_batch_jax = np.asarray(
+        jax_injection_outage_dx(
+            jacobian_inv_transposed=jnp.asarray(jacobian_inv.T),
+            outage_injection_indices=jnp.asarray(outage_batch),
+            injection_to_bus=jnp.asarray(injection_to_bus),
+            injection_active_power=jnp.asarray(injection_active_power),
+            injection_reactive_power=jnp.asarray(injection_reactive_power),
+            angle_component_indices=jnp.asarray(angle_component_indices),
+            magnitude_component_indices=jnp.asarray(magnitude_component_indices),
+        )
+    )
+    dx_batch_numpy = numpy_injection_outage_dx(
+        jacobian_inv=jacobian_inv,
+        outage_injection_indices=outage_batch,
+        injection_to_bus=injection_to_bus,
+        injection_active_power=injection_active_power,
+        injection_reactive_power=injection_reactive_power,
+        angle_component_indices=angle_component_indices,
+        magnitude_component_indices=magnitude_component_indices,
+    )
+
+    np.testing.assert_allclose(dx_batch_jax, dx_batch_numpy, rtol=TOL, atol=TOL)
+
+
+@pytest.mark.parametrize("get_net", powsybl_networks + pandapower_networks)
+def test_non_voltage_regulating_injection_outage_jax(get_net) -> None:
+    net, _, dynamic_info, _, jacobian_data = _load_test_grid(get_net)
+    injections = _get_injections_powsybl(net).reset_index(drop=True)
+    outage_p = injections["setpoint_p"].fillna(injections["p"]).to_numpy(dtype=float)
+    outage_q = injections["setpoint_q"].fillna(injections["q"]).to_numpy(dtype=float)
+
+    assert len(injections) == dynamic_info.n_injections
+
+    outage_rows = _supported_injection_outage_rows(net, injections)
+    assert outage_rows, f"No supported injection outages found for {get_net.__name__}"
+
+    outage_batch = _pad_outage_rows(outage_rows)
+    dx_batch_jax = np.asarray(
+        jax_injection_outage_dx(
+            jacobian_inv_transposed=jnp.asarray(jacobian_data.inverse_jacobian.T),
+            outage_injection_indices=jnp.asarray(outage_batch),
+            injection_to_bus=jnp.asarray(dynamic_info.injection_to_bus),
+            injection_active_power=jnp.asarray(outage_p),
+            injection_reactive_power=jnp.asarray(outage_q),
+            angle_component_indices=jnp.asarray(jacobian_data.angle_component_indices),
+            magnitude_component_indices=jnp.asarray(jacobian_data.magnitude_component_indices),
+        )
+    )
+    dx_batch_numpy = numpy_injection_outage_dx(
+        jacobian_inv=jacobian_data.inverse_jacobian,
+        outage_injection_indices=outage_batch,
+        injection_to_bus=dynamic_info.injection_to_bus,
+        injection_active_power=outage_p,
+        injection_reactive_power=outage_q,
+        angle_component_indices=jacobian_data.angle_component_indices,
+        magnitude_component_indices=jacobian_data.magnitude_component_indices,
+    )
+
+    assert dx_batch_jax.shape == (len(outage_rows), jacobian_data.inverse_jacobian.shape[0])
+    # the numpy code is tested against powsybl -> we do not need to repeat the test against powsybl
+    np.testing.assert_allclose(dx_batch_jax, dx_batch_numpy, rtol=TOL, atol=TOL)
+
+
+@pytest.mark.parametrize("get_net", powsybl_networks + pandapower_networks)
+def test_non_voltage_regulating_injection_outage_monitored_jax(get_net) -> None:
+    net, _, dynamic_info, _, jacobian_data = _load_test_grid(get_net)
+    injections = _get_injections_powsybl(net).reset_index(drop=True)
+    outage_p = injections["setpoint_p"].fillna(injections["p"]).to_numpy(dtype=float)
+    outage_q = injections["setpoint_q"].fillna(injections["q"]).to_numpy(dtype=float)
+    outage_rows = _supported_injection_outage_rows(net, injections, max_rows=5)
+    assert outage_rows, f"No supported injection outages found for {get_net.__name__}"
+    outage_batch = _pad_outage_rows(outage_rows)
+
+    branch_connected = np.asarray(dynamic_info.branch_connected, dtype=bool)
+    monitor_branch_indices = np.flatnonzero(branch_connected)[: min(5, int(branch_connected.sum()))]
+    if monitor_branch_indices.size == 0:
+        pytest.skip("No connected monitored branches available")
+
+    branch_from = np.asarray(dynamic_info.branch_from_bus, dtype=np.int64)
+    branch_to = np.asarray(dynamic_info.branch_to_bus, dtype=np.int64)
+    monitor_bus_indices = np.unique(
+        np.concatenate([branch_from[monitor_branch_indices], branch_to[monitor_branch_indices]])
+    ).astype(np.int64)
+    bus_to_mon_index = np.full(dynamic_info.n_buses, -1, dtype=np.int64)
+    bus_to_mon_index[monitor_bus_indices] = np.arange(monitor_bus_indices.size, dtype=np.int64)
+
+    monitor_theta_jax, monitor_vm_jax = non_voltage_regulating_injection_outage_monitor_buses(
+        jacobian_inv_transposed=jnp.asarray(jacobian_data.inverse_jacobian.T),
+        outage_injection_indices=jnp.asarray(outage_batch),
+        injection_to_bus=jnp.asarray(dynamic_info.injection_to_bus),
+        injection_active_power=jnp.asarray(outage_p),
+        injection_reactive_power=jnp.asarray(outage_q),
+        angle_component_indices=jnp.asarray(jacobian_data.angle_component_indices),
+        magnitude_component_indices=jnp.asarray(jacobian_data.magnitude_component_indices),
+        monitor_bus_indices=jnp.asarray(monitor_bus_indices),
+        v_mag_hat=jnp.asarray(dynamic_info.bus_voltage_magnitudes),
+        theta_hat=jnp.asarray(dynamic_info.bus_voltage_angles_rad),
+    )
+    monitored_results = non_voltage_regulating_injection_outage_monitored(
+        jacobian_inv_transposed=jnp.asarray(jacobian_data.inverse_jacobian.T),
+        outage_injection_indices=jnp.asarray(outage_batch),
+        injection_to_bus=jnp.asarray(dynamic_info.injection_to_bus),
+        injection_active_power=jnp.asarray(outage_p),
+        injection_reactive_power=jnp.asarray(outage_q),
+        angle_component_indices=jnp.asarray(jacobian_data.angle_component_indices),
+        magnitude_component_indices=jnp.asarray(jacobian_data.magnitude_component_indices),
+        monitor_bus_indices=jnp.asarray(monitor_bus_indices),
+        v_mag_hat=jnp.asarray(dynamic_info.bus_voltage_magnitudes),
+        theta_hat=jnp.asarray(dynamic_info.bus_voltage_angles_rad),
+        branch_from=jnp.asarray(branch_from),
+        branch_to=jnp.asarray(branch_to),
+        y_ff=jnp.asarray(dynamic_info.branch_effective_admittance_from_from),
+        y_ft=jnp.asarray(dynamic_info.branch_effective_admittance_from_to),
+        y_tf=jnp.asarray(dynamic_info.branch_effective_admittance_to_from),
+        y_tt=jnp.asarray(dynamic_info.branch_effective_admittance_to_to),
+        monitor_branch_indices=jnp.asarray(monitor_branch_indices),
+        bus_to_mon_index=jnp.asarray(bus_to_mon_index),
+    )
+    dx_batch_numpy = numpy_injection_outage_dx(
+        jacobian_inv=jacobian_data.inverse_jacobian,
+        outage_injection_indices=outage_batch,
+        injection_to_bus=dynamic_info.injection_to_bus,
+        injection_active_power=outage_p,
+        injection_reactive_power=outage_q,
+        angle_component_indices=jacobian_data.angle_component_indices,
+        magnitude_component_indices=jacobian_data.magnitude_component_indices,
+    )
+
+    expected_theta = np.broadcast_to(
+        np.asarray(dynamic_info.bus_voltage_angles_rad, dtype=float)[monitor_bus_indices],
+        (len(outage_rows), monitor_bus_indices.size),
+    ).copy()
+    expected_vm = np.broadcast_to(
+        np.asarray(dynamic_info.bus_voltage_magnitudes, dtype=float)[monitor_bus_indices],
+        (len(outage_rows), monitor_bus_indices.size),
+    ).copy()
+    theta_component_indices = np.asarray(jacobian_data.angle_component_indices, dtype=np.int64)[monitor_bus_indices]
+    vm_component_indices = np.asarray(jacobian_data.magnitude_component_indices, dtype=np.int64)[monitor_bus_indices]
+    theta_mask = theta_component_indices >= 0
+    vm_mask = vm_component_indices >= 0
+    expected_theta[:, theta_mask] += dx_batch_numpy[:, theta_component_indices[theta_mask]]
+    expected_vm[:, vm_mask] += dx_batch_numpy[:, vm_component_indices[vm_mask]]
+
+    np.testing.assert_allclose(np.asarray(monitor_theta_jax), expected_theta, rtol=TOL, atol=TOL)
+    np.testing.assert_allclose(np.asarray(monitor_vm_jax), expected_vm, rtol=TOL, atol=TOL)
+    np.testing.assert_allclose(np.asarray(monitored_results.n_1_theta), expected_theta, rtol=TOL, atol=TOL)
+    np.testing.assert_allclose(np.asarray(monitored_results.n_1_voltage), expected_vm, rtol=TOL, atol=TOL)
+
+    f_pos = bus_to_mon_index[branch_from[monitor_branch_indices]]
+    t_pos = bus_to_mon_index[branch_to[monitor_branch_indices]]
+    assert np.all(f_pos >= 0)
+    assert np.all(t_pos >= 0)
+
+    v_post = expected_vm * np.exp(1j * expected_theta)
+    y_ff_mon = np.asarray(dynamic_info.branch_effective_admittance_from_from)[monitor_branch_indices]
+    y_ft_mon = np.asarray(dynamic_info.branch_effective_admittance_from_to)[monitor_branch_indices]
+    y_tf_mon = np.asarray(dynamic_info.branch_effective_admittance_to_from)[monitor_branch_indices]
+    y_tt_mon = np.asarray(dynamic_info.branch_effective_admittance_to_to)[monitor_branch_indices]
+    v_from = v_post[:, f_pos]
+    v_to = v_post[:, t_pos]
+    i_from = y_ff_mon[None, :] * v_from + y_ft_mon[None, :] * v_to
+    i_to = y_tf_mon[None, :] * v_from + y_tt_mon[None, :] * v_to
+    s_from = v_from * np.conj(i_from)
+    s_to = v_to * np.conj(i_to)
+
+    np.testing.assert_allclose(np.asarray(monitored_results.n_1_i_from), i_from, rtol=1e-9, atol=1e-9)
+    np.testing.assert_allclose(np.asarray(monitored_results.n_1_i_to), i_to, rtol=1e-9, atol=1e-9)
+    np.testing.assert_allclose(np.asarray(monitored_results.n_1_p_from), s_from.real, rtol=1e-9, atol=1e-9)
+    np.testing.assert_allclose(np.asarray(monitored_results.n_1_p_to), s_to.real, rtol=1e-9, atol=1e-9)
+    np.testing.assert_allclose(np.asarray(monitored_results.n_1_q_from), s_from.imag, rtol=1e-9, atol=1e-9)
+    np.testing.assert_allclose(np.asarray(monitored_results.n_1_q_to), s_to.imag, rtol=1e-9, atol=1e-9)
diff --git a/tests/numpy/test_injection_outage.py b/tests/numpy/test_injection_outage.py
new file mode 100644
index 0000000..054bb21
--- /dev/null
+++ b/tests/numpy/test_injection_outage.py
@@ -0,0 +1,164 @@
+# Copyright 2026 50Hertz Transmission GmbH and Elia Transmission Belgium SA/NV
+#
+# This Source Code Form is subject to the terms of the Mozilla Public License, v. 2.0.
+# If a copy of the MPL was not distributed with this file,
+# you can obtain one at https://mozilla.org/MPL/2.0/.
+# Mozilla Public License, version 2.0
+
+from copy import deepcopy
+
+import numpy as np
+import pypowsybl
+import pytest
+
+from dc_plus.example_grids.pypowsbl.example_grids import POWSYBL_NETWORKS_SHORT_LIST
+from dc_plus.importing.powsybl.powsybl_import import DANGLING_BUS_STRING_SUFFIX, _get_injections_powsybl
+from dc_plus.importing.powsybl.powsybl_loadflow_parameter import get_powsybl_loadflow_parameter
+from dc_plus.importing.powsybl.powsybl_network_helpers import _load_test_grid
+from dc_plus.numpy.injection_outage import non_voltage_regulating_injection_outage_dx
+from dc_plus.preprocess.create_network_data import create_network_data
+
+powsybl_networks = POWSYBL_NETWORKS_SHORT_LIST
+TOL = 1e-10
+
+
+def _connected_injection_indices(injections):
+    candidates = np.flatnonzero(injections["connected"].to_numpy(dtype=bool))
+    if candidates.size == 0:
+        pytest.skip("No connected injections available for an injection outage test")
+    return candidates
+
+
+def _expected_bus_delta(dynamic_info, outage_injection_indices, outage_p, outage_q):
+    delta_p = np.zeros(dynamic_info.n_buses, dtype=float)
+    delta_q = np.zeros(dynamic_info.n_buses, dtype=float)
+    outage_indices = np.asarray(outage_injection_indices, dtype=np.int64)
+    valid_mask = outage_indices >= 0
+    if not np.any(valid_mask):
+        return delta_p, delta_q
+
+    valid_outages = outage_indices[valid_mask]
+    buses = dynamic_info.injection_to_bus[valid_outages]
+    np.add.at(delta_p, buses, -outage_p[valid_outages])
+    np.add.at(delta_q, buses, -outage_q[valid_outages])
+    return delta_p, delta_q
+
+
+def _pad_outage_rows(outage_rows: list[np.ndarray]) -> np.ndarray:
+    if not outage_rows:
+        return np.empty((0, 0), dtype=np.int64)
+
+    max_outages = max(outage_row.size for outage_row in outage_rows)
+    padded = np.full((len(outage_rows), max_outages), -1, dtype=np.int64)
+    for row_idx, outage_row in enumerate(outage_rows):
+        padded[row_idx, : outage_row.size] = outage_row
+    return padded
+
+
+@pytest.mark.parametrize("get_net", powsybl_networks)
+def test_non_voltage_regulating_injection_outage_numpy(get_net) -> None:
+    """Test non voltage regulating injection outages with a one-step update against powsybl load flow results.
+
+    This does not test HVDC converter and dangling line outages.
+    It does test voltage regulating generator outages, as the powsybl one step does also not switch from PV to PQ
+    node. This therefore does not account for the voltage drop of the disconnected generator.
+
+    In DC+ dangling lines are modelled as an additional bus with a fixed injection p0 and q0
+    a one sided line is added between the original bus and the new bus to account for the line losses
+    of the dangling line. Here numerical errors occur.
+    The test still verifies the voltages of the bus where the dangling line is attached to with the
+    full precision.
+
+    Parameters
+    ----------
+    get_net : Callable[[], pypowsybl.network.Network]
+        A fixture that returns a powsybl network to test on.
+    """
+    net, _, dynamic_info, string_info, jacobian_data = _load_test_grid(get_net)
+    injections = _get_injections_powsybl(net).reset_index(drop=True)
+    injection_powsybl_type = net.get_injections(attributes=["type"])
+    outage_p = injections["setpoint_p"].fillna(injections["p"]).to_numpy(dtype=float)
+    outage_q = injections["setpoint_q"].fillna(injections["q"]).to_numpy(dtype=float)
+    dangling_bus_mask = np.char.endswith(string_info.bus_ids.astype(str), DANGLING_BUS_STRING_SUFFIX)
+
+    assert len(injections) == dynamic_info.n_injections
+
+    accepted_cases: list[tuple[str, np.ndarray, object]] = []
+    for injection_idx in _connected_injection_indices(injections):
+        injection_id = injections.loc[injection_idx, "id_str"]
+        if (injection_powsybl_type.loc[injection_id, "type"]) in (["HVDC_CONVERTER_STATION", "DANGLING_LINE"]):
+            # - HVDC would need a multi outage
+            # - Dangling line removal would change the network for DC+
+            continue
+
+        net_n1 = deepcopy(net)
+        net_n1.remove_elements(injection_id)
+
+        loadflow_parameter = get_powsybl_loadflow_parameter("one_step")
+        loadflow_res = pypowsybl.loadflow.run_ac(net_n1, parameters=loadflow_parameter)[0]
+        assert loadflow_res.iteration_count <= 1
+        if loadflow_res.status != pypowsybl._pypowsybl.LoadFlowComponentStatus.CONVERGED:
+            continue
+
+        _, dynamic_info_n1, string_info_n1 = create_network_data(net_n1)
+        if not np.array_equal(string_info.bus_ids, string_info_n1.bus_ids):
+            continue
+
+        same_angle_structure = np.array_equal(
+            dynamic_info.pvpq_buses_indices_pvpq_order,
+            dynamic_info_n1.pvpq_buses_indices_pvpq_order,
+        )
+        same_magnitude_structure = np.array_equal(
+            dynamic_info.pq_buses_indices,
+            dynamic_info_n1.pq_buses_indices,
+        )
+        if not (same_angle_structure and same_magnitude_structure):
+            continue
+
+        outage_row = np.array([injection_idx], dtype=np.int64)
+        expected_delta_p, expected_delta_q = _expected_bus_delta(dynamic_info, outage_row, outage_p, outage_q)
+        actual_delta_p = dynamic_info_n1.bus_active_power - dynamic_info.bus_active_power
+        actual_delta_q = dynamic_info_n1.bus_reactive_power - dynamic_info.bus_reactive_power
+        p_buses = dynamic_info.pvpq_buses_indices_pvpq_order
+        q_buses = dynamic_info.pq_buses_indices
+        if not (
+            np.allclose(actual_delta_p[p_buses], expected_delta_p[p_buses], rtol=TOL, atol=TOL)
+            and np.allclose(actual_delta_q[q_buses], expected_delta_q[q_buses], rtol=TOL, atol=TOL)
+        ):
+            continue
+
+        accepted_cases.append((injection_id, outage_row, dynamic_info_n1))
+
+    assert accepted_cases, f"No injection outages matched the fixed-injection assumptions for {get_net.__name__}"
+
+    outage_batch = _pad_outage_rows([outage_row for _, outage_row, _ in accepted_cases])
+    dx_batch = non_voltage_regulating_injection_outage_dx(
+        jacobian_inv=jacobian_data.inverse_jacobian,
+        outage_injection_indices=outage_batch,
+        injection_to_bus=dynamic_info.injection_to_bus,
+        injection_active_power=outage_p,
+        injection_reactive_power=outage_q,
+        angle_component_indices=jacobian_data.angle_component_indices,
+        magnitude_component_indices=jacobian_data.magnitude_component_indices,
+    )
+
+    assert dx_batch.shape == (len(accepted_cases), jacobian_data.inverse_jacobian.shape[0])
+
+    for case_idx, (_, _, dynamic_info_n1) in enumerate(accepted_cases):
+        dx = dx_batch[case_idx]
+        theta_updated = dynamic_info.bus_voltage_angles_rad.copy()
+        vm_updated = dynamic_info.bus_voltage_magnitudes.copy()
+        theta_updated[dynamic_info.pvpq_buses_indices_pvpq_order] += dx[jacobian_data.is_angle_component]
+        vm_updated[dynamic_info.pq_buses_indices] += dx[jacobian_data.is_magnitude_component]
+        np.testing.assert_allclose(
+            theta_updated[~dangling_bus_mask],
+            dynamic_info_n1.bus_voltage_angles_rad[~dangling_bus_mask],
+            rtol=TOL,
+            atol=TOL,
+        )
+        np.testing.assert_allclose(
+            vm_updated[~dangling_bus_mask],
+            dynamic_info_n1.bus_voltage_magnitudes[~dangling_bus_mask],
+            rtol=TOL,
+            atol=TOL,
+        )