diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/main_core/math_output_queries.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/main_core/math_output_queries.hpp
index 15d30677d6..931f67a3d3 100644
--- a/power_grid_model_c/power_grid_model/include/power_grid_model/main_core/math_output_queries.hpp
+++ b/power_grid_model_c/power_grid_model/include/power_grid_model/main_core/math_output_queries.hpp
@@ -19,6 +19,29 @@
 
 namespace power_grid_model::main_core {
 
+template <solver_output_type SolverOutputType>
+constexpr auto get_voltage_output(MathOutput<std::vector<SolverOutputType>> const& math_output, Idx2D const& math_id) {
+    if constexpr (steady_state_solver_output_type<SolverOutputType>) {
+        return math_output.solver_output[math_id.group].u[math_id.pos];
+    } else if constexpr (short_circuit_solver_output_type<SolverOutputType>) {
+        return math_output.solver_output[math_id.group].u_bus[math_id.pos];
+    } else {
+        static_assert(false, "Unsupported solver output type for voltage output retrieval");
+    }
+}
+
+template <steady_state_solver_output_type SolverOutputType>
+constexpr auto get_bus_injection_output_from_topo_id(MathOutput<std::vector<SolverOutputType>> const& math_output,
+                                                     Idx2D const& topo_id) {
+    return math_output.supernode_output[topo_id.group].bus_injection[topo_id.pos];
+}
+
+template <steady_state_solver_output_type SolverOutputType>
+constexpr auto get_bus_injection_output_from_math_id(MathOutput<std::vector<SolverOutputType>> const& math_output,
+                                                     Idx2D const& math_id) {
+    return math_output.solver_output[math_id.group].bus_injection[math_id.pos];
+}
+
 template <typename Component, solver_output_type SolverOutputType>
 constexpr auto const& get_component_output(MathOutput<std::vector<SolverOutputType>> const& math_output,
                                            Idx2D math_id) {
diff --git a/power_grid_model_c/power_grid_model/include/power_grid_model/main_core/output.hpp b/power_grid_model_c/power_grid_model/include/power_grid_model/main_core/output.hpp
index 781fb84eb4..e008e8b733 100644
--- a/power_grid_model_c/power_grid_model/include/power_grid_model/main_core/output.hpp
+++ b/power_grid_model_c/power_grid_model/include/power_grid_model/main_core/output.hpp
@@ -5,6 +5,7 @@
 #pragma once
 
 #include "container_queries.hpp"
+#include "math_output_queries.hpp"
 #include "state.hpp"
 #include "state_queries.hpp"
 
@@ -111,8 +112,8 @@ constexpr auto output_result(Component const& node, MainModelState<ComponentCont
         return node.template get_null_output<sym>();
     }
 
-    return node.template get_output<sym>(math_output.solver_output[math_id.group].u[math_id.pos],
-                                         math_output.supernode_output[topo_id.group].bus_injection[topo_id.pos]);
+    return node.template get_output<sym>(get_voltage_output(math_output, math_id),
+                                         get_bus_injection_output_from_topo_id(math_output, topo_id));
 }
 template <std::derived_from<Node> Component, class ComponentContainer,
           short_circuit_solver_output_type SolverOutputType>
@@ -125,7 +126,7 @@ constexpr auto output_result(Component const& node, MainModelState<ComponentCont
         return node.get_null_sc_output();
     }
 
-    return node.get_sc_output(math_output.solver_output[math_id.group].u_bus[math_id.pos]);
+    return node.get_sc_output(get_voltage_output(math_output, math_id));
 }
 
 // output link
@@ -159,28 +160,29 @@ inline auto output_result(Component const& link, MainModelState<ComponentContain
 // output branch
 template <std::derived_from<Edge> Component, steady_state_solver_output_type SolverOutputType>
     requires(!std::same_as<Component, Link>) // TODO(mgovers): cleanup v2: change back to only derived_from<Branch>
-constexpr auto output_result(Component const& branch, std::vector<SolverOutputType> const& solver_output,
+constexpr auto output_result(Component const& branch, MathOutput<std::vector<SolverOutputType>> const& math_output,
                              Idx2D math_id) {
     using sym = decode_symmetry_v<SolverOutputType>;
 
     if (math_id.group == disconnected) {
         return branch.template get_null_output<sym>();
     }
-    return branch.template get_output<sym>(solver_output[math_id.group].branch[math_id.pos]);
+    return branch.template get_output<sym>(get_component_output<Component>(math_output, math_id));
 }
 // TODO(mgovers): cleanup v2: change back to only derived_from<Branch>
 template <std::derived_from<Edge> Component, short_circuit_solver_output_type SolverOutputType>
     requires(!std::same_as<Component, Link>) // TODO(mgovers): cleanup v2: change back to only derived_from<Branch>
-inline auto output_result(Component const& branch, std::vector<SolverOutputType> const& solver_output, Idx2D math_id) {
+inline auto output_result(Component const& branch, MathOutput<std::vector<SolverOutputType>> const& math_output,
+                          Idx2D const& math_id) {
     if (math_id.group == disconnected) {
         return branch.get_null_sc_output();
     }
-    return branch.get_sc_output(solver_output[math_id.group].branch[math_id.pos]);
+    return branch.get_sc_output(get_component_output<Component>(math_output, math_id));
 }
 
 // output branch3
 template <std::derived_from<Branch3> Component, steady_state_solver_output_type SolverOutputType>
-constexpr auto output_result(Component const& branch3, std::vector<SolverOutputType> const& solver_output,
+constexpr auto output_result(Component const& branch3, MathOutput<std::vector<SolverOutputType>> const& math_output,
                              Idx2DBranch3 const& math_id) {
     using sym = decode_symmetry_v<SolverOutputType>;
 
@@ -188,76 +190,79 @@ constexpr auto output_result(Component const& branch3, std::vector<SolverOutputT
         return branch3.template get_null_output<sym>();
     }
 
-    auto const& branches = solver_output[math_id.group].branch;
-    return branch3.template get_output<sym>(branches[math_id.pos[0]], branches[math_id.pos[1]],
-                                            branches[math_id.pos[2]]);
+    return branch3.template get_output<sym>(
+        get_component_output<Branch3>(math_output, {.group = math_id.group, .pos = math_id.pos[0]}),
+        get_component_output<Branch3>(math_output, {.group = math_id.group, .pos = math_id.pos[1]}),
+        get_component_output<Branch3>(math_output, {.group = math_id.group, .pos = math_id.pos[2]}));
 }
 template <std::derived_from<Branch3> Component, short_circuit_solver_output_type SolverOutputType>
-inline auto output_result(Component const& branch3, std::vector<SolverOutputType> const& solver_output,
+inline auto output_result(Component const& branch3, MathOutput<std::vector<SolverOutputType>> const& math_output,
                           Idx2DBranch3 const& math_id) {
     if (math_id.group == disconnected) {
         return branch3.get_null_sc_output();
     }
 
-    auto const& branches = solver_output[math_id.group].branch;
-    return branch3.get_sc_output(branches[math_id.pos[0]], branches[math_id.pos[1]], branches[math_id.pos[2]]);
+    return branch3.get_sc_output(
+        get_component_output<Branch3>(math_output, {.group = math_id.group, .pos = math_id.pos[0]}),
+        get_component_output<Branch3>(math_output, {.group = math_id.group, .pos = math_id.pos[1]}),
+        get_component_output<Branch3>(math_output, {.group = math_id.group, .pos = math_id.pos[2]}));
 }
 
 // output source
 template <std::derived_from<Source> Component, steady_state_solver_output_type SolverOutputType>
-constexpr auto output_result(Component const& source, std::vector<SolverOutputType> const& solver_output,
+constexpr auto output_result(Component const& source, MathOutput<std::vector<SolverOutputType>> const& math_output,
                              Idx2D const& math_id) {
     using sym = decode_symmetry_v<SolverOutputType>;
 
     if (math_id.group == disconnected) {
         return source.template get_null_output<sym>();
     }
-    return source.template get_output<sym>(solver_output[math_id.group].source[math_id.pos]);
+    return source.template get_output<sym>(get_component_output<Source>(math_output, math_id));
 }
 template <std::derived_from<Source> Component, short_circuit_solver_output_type SolverOutputType>
-inline auto output_result(Component const& source, std::vector<SolverOutputType> const& solver_output,
+inline auto output_result(Component const& source, MathOutput<std::vector<SolverOutputType>> const& math_output,
                           Idx2D const& math_id) {
     if (math_id.group == disconnected) {
         return source.get_null_sc_output();
     }
-    return source.get_sc_output(solver_output[math_id.group].source[math_id.pos]);
+    return source.get_sc_output(get_component_output<Source>(math_output, math_id));
 }
 
 // output load gen
 template <std::derived_from<GenericLoadGen> Component, steady_state_solver_output_type SolverOutputType>
-constexpr auto output_result(Component const& load_gen, std::vector<SolverOutputType> const& solver_output,
+constexpr auto output_result(Component const& load_gen, MathOutput<std::vector<SolverOutputType>> const& math_output,
                              Idx2D const& math_id) {
     using sym = decode_symmetry_v<SolverOutputType>;
 
     if (math_id.group == disconnected) {
         return load_gen.template get_null_output<sym>();
     }
-    return load_gen.template get_output<sym>(solver_output[math_id.group].load_gen[math_id.pos]);
+    return load_gen.template get_output<sym>(get_component_output<GenericLoadGen>(math_output, math_id));
 }
 template <std::derived_from<GenericLoadGen> Component, short_circuit_solver_output_type SolverOutputType>
-inline auto output_result(Component const& load_gen, std::vector<SolverOutputType> const& /*solver_output*/,
+inline auto output_result(Component const& load_gen, MathOutput<std::vector<SolverOutputType>> const& /*math_output*/,
                           Idx2D const& /*math_id*/) {
     return load_gen.get_null_sc_output();
 }
 
 // output shunt
 template <std::derived_from<Shunt> Component, steady_state_solver_output_type SolverOutputType>
-constexpr auto output_result(Component const& shunt, std::vector<SolverOutputType> const& solver_output,
+constexpr auto output_result(Component const& shunt, MathOutput<std::vector<SolverOutputType>> const& math_output,
                              Idx2D const& math_id) {
     using sym = decode_symmetry_v<SolverOutputType>;
 
     if (math_id.group == disconnected) {
         return shunt.template get_null_output<sym>();
     }
-    return shunt.template get_output<sym>(solver_output[math_id.group].shunt[math_id.pos]);
+    return shunt.template get_output<sym>(get_component_output<Shunt>(math_output, math_id));
 }
 template <std::derived_from<Shunt> Component, short_circuit_solver_output_type SolverOutputType>
-inline auto output_result(Component const& shunt, std::vector<SolverOutputType> const& solver_output,
+inline auto output_result(Component const& shunt, MathOutput<std::vector<SolverOutputType>> const& math_output,
                           Idx2D const& math_id) {
     if (math_id.group == disconnected) {
         return shunt.get_null_sc_output();
     }
-    return shunt.get_sc_output(solver_output[math_id.group].shunt[math_id.pos]);
+    return shunt.get_sc_output(get_component_output<Shunt>(math_output, math_id));
 }
 
 // output voltage sensor
@@ -265,20 +270,21 @@ template <std::derived_from<GenericVoltageSensor> Component, class ComponentCont
           steady_state_solver_output_type SolverOutputType>
     requires model_component_state_c<MainModelState, ComponentContainer, Component>
 constexpr auto output_result(Component const& voltage_sensor, MainModelState<ComponentContainer> const& state,
-                             std::vector<SolverOutputType> const& solver_output, Idx const node_seq) {
+                             MathOutput<std::vector<SolverOutputType>> const& math_output, Idx const node_seq) {
     using sym = decode_symmetry_v<SolverOutputType>;
 
     Idx2D const node_math_id = state.topo_comp_coup->node[node_seq];
     if (node_math_id.group == disconnected) {
         return voltage_sensor.template get_null_output<sym>();
     }
-    return voltage_sensor.template get_output<sym>(solver_output[node_math_id.group].u[node_math_id.pos]);
+    return voltage_sensor.template get_output<sym>(get_voltage_output(math_output, node_math_id));
 }
 template <std::derived_from<GenericVoltageSensor> Component, class ComponentContainer,
           short_circuit_solver_output_type SolverOutputType>
     requires model_component_state_c<MainModelState, ComponentContainer, Component>
 inline auto output_result(Component const& voltage_sensor, MainModelState<ComponentContainer> const& /* state */,
-                          std::vector<SolverOutputType> const& /* solver_output */, Idx const /* node_seq */) {
+                          MathOutput<std::vector<SolverOutputType>> const& /* math_output */,
+                          Idx const /* node_seq */) {
     return voltage_sensor.get_null_sc_output();
 }
 
@@ -287,7 +293,7 @@ template <std::derived_from<GenericPowerSensor> Component, class ComponentContai
           steady_state_solver_output_type SolverOutputType>
     requires model_component_state_c<MainModelState, ComponentContainer, Component>
 constexpr auto output_result(Component const& power_sensor, MainModelState<ComponentContainer> const& state,
-                             std::vector<SolverOutputType> const& solver_output, Idx const obj_seq) {
+                             MathOutput<std::vector<SolverOutputType>> const& math_output, Idx const obj_seq) {
     using sym = decode_symmetry_v<SolverOutputType>;
 
     auto const terminal_type = power_sensor.get_terminal_type();
@@ -336,19 +342,19 @@ constexpr auto output_result(Component const& power_sensor, MainModelState<Compo
     case branch3_2:
         [[fallthrough]];
     case branch3_3:
-        return power_sensor.template get_output<sym>(solver_output[obj_math_id.group].branch[obj_math_id.pos].s_f);
+        return power_sensor.template get_output<sym>(get_component_output<Branch>(math_output, obj_math_id).s_f);
     case branch_to:
-        return power_sensor.template get_output<sym>(solver_output[obj_math_id.group].branch[obj_math_id.pos].s_t);
+        return power_sensor.template get_output<sym>(get_component_output<Branch>(math_output, obj_math_id).s_t);
     case source:
-        return power_sensor.template get_output<sym>(solver_output[obj_math_id.group].source[obj_math_id.pos].s);
+        return power_sensor.template get_output<sym>(get_component_output<Source>(math_output, obj_math_id).s);
     case shunt:
-        return power_sensor.template get_output<sym>(solver_output[obj_math_id.group].shunt[obj_math_id.pos].s);
+        return power_sensor.template get_output<sym>(get_component_output<Shunt>(math_output, obj_math_id).s);
     case load:
         [[fallthrough]];
     case generator:
-        return power_sensor.template get_output<sym>(solver_output[obj_math_id.group].load_gen[obj_math_id.pos].s);
+        return power_sensor.template get_output<sym>(get_component_output<GenericLoadGen>(math_output, obj_math_id).s);
     case node:
-        return power_sensor.template get_output<sym>(solver_output[obj_math_id.group].bus_injection[obj_math_id.pos]);
+        return power_sensor.template get_output<sym>(get_bus_injection_output_from_math_id(math_output, obj_math_id));
     default:
         throw MissingCaseForEnumError{std::format("{} output_result()", Component::name), terminal_type};
     }
@@ -356,9 +362,9 @@ constexpr auto output_result(Component const& power_sensor, MainModelState<Compo
 template <std::derived_from<GenericPowerSensor> Component, class ComponentContainer,
           short_circuit_solver_output_type SolverOutputType>
     requires model_component_state_c<MainModelState, ComponentContainer, Component>
-constexpr auto output_result(Component const& power_or_current_sensor,
-                             MainModelState<ComponentContainer> const& /* state */,
-                             std::vector<SolverOutputType> const& /* solver_output */, Idx const /* obj_seq */) {
+constexpr auto
+output_result(Component const& power_or_current_sensor, MainModelState<ComponentContainer> const& /* state */,
+              MathOutput<std::vector<SolverOutputType>> const& /* math_output */, Idx const /* obj_seq */) {
     return power_or_current_sensor.get_null_sc_output();
 }
 
@@ -367,7 +373,7 @@ template <std::derived_from<GenericCurrentSensor> Component, class ComponentCont
           steady_state_solver_output_type SolverOutputType>
     requires model_component_state_c<MainModelState, ComponentContainer, Component>
 constexpr auto output_result(Component const& current_sensor, MainModelState<ComponentContainer> const& state,
-                             std::vector<SolverOutputType> const& solver_output, Idx const obj_seq) {
+                             MathOutput<std::vector<SolverOutputType>> const& math_output, Idx const obj_seq) {
     using sym = decode_symmetry_v<SolverOutputType>;
 
     auto const terminal_type = current_sensor.get_terminal_type();
@@ -411,11 +417,11 @@ constexpr auto output_result(Component const& current_sensor, MainModelState<Com
     case branch3_2:
         [[fallthrough]];
     case branch3_3:
-        return current_sensor.template get_output<sym>(solver_output[obj_math_id.group].branch[obj_math_id.pos].i_f,
-                                                       solver_output[node_from_math_id.group].u[node_from_math_id.pos]);
+        return current_sensor.template get_output<sym>(get_component_output<Branch>(math_output, obj_math_id).i_f,
+                                                       get_voltage_output(math_output, node_from_math_id));
     case branch_to:
-        return current_sensor.template get_output<sym>(solver_output[obj_math_id.group].branch[obj_math_id.pos].i_t,
-                                                       solver_output[node_to_math_id.group].u[node_to_math_id.pos]);
+        return current_sensor.template get_output<sym>(get_component_output<Branch>(math_output, obj_math_id).i_t,
+                                                       get_voltage_output(math_output, node_to_math_id));
     default:
         throw MissingCaseForEnumError{std::format("{} output_result()", Component::name), terminal_type};
     }
@@ -423,9 +429,9 @@ constexpr auto output_result(Component const& current_sensor, MainModelState<Com
 template <std::derived_from<GenericCurrentSensor> Component, class ComponentContainer,
           short_circuit_solver_output_type SolverOutputType>
     requires model_component_state_c<MainModelState, ComponentContainer, Component>
-constexpr auto output_result(Component const& power_or_current_sensor,
-                             MainModelState<ComponentContainer> const& /* state */,
-                             std::vector<SolverOutputType> const& /* solver_output */, Idx const /* obj_seq */) {
+constexpr auto
+output_result(Component const& power_or_current_sensor, MainModelState<ComponentContainer> const& /* state */,
+              MathOutput<std::vector<SolverOutputType>> const& /* math_output */, Idx const /* obj_seq */) {
     return power_or_current_sensor.get_null_sc_output();
 }
 
@@ -435,7 +441,7 @@ template <std::derived_from<Fault> Component, class ComponentContainer,
     requires model_component_state_c<MainModelState, ComponentContainer, Component> &&
              model_component_state_c<MainModelState, ComponentContainer, Node>
 constexpr auto output_result(Component const& fault, MainModelState<ComponentContainer> const& /* state */,
-                             std::vector<SolverOutputType> const& /* solver_output */, Idx2D /* math_id */) {
+                             MathOutput<std::vector<SolverOutputType>> const& /* math_output */, Idx2D /* math_id */) {
     return fault.get_output();
 }
 template <std::derived_from<Fault> Component, class ComponentContainer,
@@ -443,13 +449,13 @@ template <std::derived_from<Fault> Component, class ComponentContainer,
     requires model_component_state_c<MainModelState, ComponentContainer, Component> &&
              model_component_state_c<MainModelState, ComponentContainer, Node>
 inline auto output_result(Component const& fault, MainModelState<ComponentContainer> const& state,
-                          std::vector<SolverOutputType> const& solver_output, Idx2D math_id) {
+                          MathOutput<std::vector<SolverOutputType>> const& math_output, Idx2D math_id) {
     if (math_id.group == disconnected) {
         return fault.get_null_sc_output();
     }
 
     auto const u_rated = get_component<Node>(state.components, fault.get_fault_object()).u_rated();
-    return fault.get_sc_output(solver_output[math_id.group].fault[math_id.pos], u_rated);
+    return fault.get_sc_output(get_component_output<Fault>(math_output, math_id), u_rated);
 }
 
 // output transformer tap regulator
@@ -519,10 +525,10 @@ constexpr void output_result(MainModelState<ComponentContainer> const& state,
                 return output_result<Component, ComponentContainer>(link, state, math_output, topo_id);
             });
     } else {
-        detail::produce_output<Component, Idx2D>(
-            state, output, [&math_output](Component const& link, Idx2D const& math_id) {
-                return output_result<Edge>(link, math_output.solver_output, math_id);
-            });
+        detail::produce_output<Component, Idx2D>(state, output,
+                                                 [&math_output](Component const& link, Idx2D const& math_id) {
+                                                     return output_result<Edge>(link, math_output, math_id);
+                                                 });
     }
 }
 
@@ -530,63 +536,66 @@ constexpr void output_result(MainModelState<ComponentContainer> const& state,
 template <std::derived_from<Base> Component, class ComponentContainer, solver_output_type SolverOutputType,
           non_owning_view_c ComponentOutput>
     requires model_component_state_c<MainModelState, ComponentContainer, Component> &&
-             requires(Component const& component, std::vector<SolverOutputType> const& solver_output, Idx2D math_id) {
+             requires(Component const& component, MathOutput<std::vector<SolverOutputType>> const& math_output,
+                      Idx2D math_id) {
                  {
-                     output_result<Component>(component, solver_output, math_id)
+                     output_result<Component>(component, math_output, math_id)
                  } -> detail::assignable_to<std::ranges::range_reference_t<ComponentOutput>>;
              }
 constexpr void output_result(MainModelState<ComponentContainer> const& state,
                              MathOutput<std::vector<SolverOutputType>> const& math_output, ComponentOutput output) {
     detail::produce_output<Component, Idx2D>(state, output, [&math_output](Component const& component, Idx2D math_id) {
-        return output_result<Component>(component, math_output.solver_output, math_id);
+        return output_result<Component>(component, math_output, math_id);
     });
 }
 template <std::derived_from<Base> Component, class ComponentContainer, solver_output_type SolverOutputType,
           non_owning_view_c ComponentOutput>
     requires model_component_state_c<MainModelState, ComponentContainer, Component> &&
              requires(Component const& component, MainModelState<ComponentContainer> const& state,
-                      std::vector<SolverOutputType> const& solver_output, Idx2D math_id) {
+                      MathOutput<std::vector<SolverOutputType>> const& math_output, Idx2D math_id) {
                  {
-                     output_result<Component>(component, state, solver_output, math_id)
+                     output_result<Component>(component, state, math_output, math_id)
                  } -> detail::assignable_to<std::ranges::range_reference_t<ComponentOutput>>;
-             }
+             } &&
+             (!std::same_as<Component, Link>) // TODO(mgovers): cleanup v2: this requirement should no longer be needed
+                                              // after link output is cleaned up
 constexpr void output_result(MainModelState<ComponentContainer> const& state,
                              MathOutput<std::vector<SolverOutputType>> const& math_output, ComponentOutput output) {
     detail::produce_output<Component, Idx2D>(
         state, output, [&state, &math_output](Component const& component, Idx2D const math_id) {
-            return output_result<Component>(component, state, math_output.solver_output, math_id);
+            return output_result<Component>(component, state, math_output, math_id);
         });
 }
 template <std::derived_from<Base> Component, class ComponentContainer, solver_output_type SolverOutputType,
           non_owning_view_c ComponentOutput>
     requires model_component_state_c<MainModelState, ComponentContainer, Component> &&
              requires(Component const& component, MainModelState<ComponentContainer> const& state,
-                      std::vector<SolverOutputType> const& solver_output, Idx obj_seq) {
+                      MathOutput<std::vector<SolverOutputType>> const& math_output, Idx obj_seq) {
                  {
-                     output_result<Component>(component, state, solver_output, obj_seq)
+                     output_result<Component>(component, state, math_output, obj_seq)
                  } -> detail::assignable_to<std::ranges::range_reference_t<ComponentOutput>>;
              }
 constexpr void output_result(MainModelState<ComponentContainer> const& state,
                              MathOutput<std::vector<SolverOutputType>> const& math_output, ComponentOutput output) {
     detail::produce_output<Component, Idx>(
         state, output, [&state, &math_output](Component const& component, Idx const obj_seq) {
-            return output_result<Component, ComponentContainer>(component, state, math_output.solver_output, obj_seq);
+            return output_result<Component, ComponentContainer>(component, state, math_output, obj_seq);
         });
 }
 template <std::derived_from<Base> Component, class ComponentContainer, solver_output_type SolverOutputType,
           non_owning_view_c ComponentOutput>
     requires model_component_state_c<MainModelState, ComponentContainer, Component> &&
-             requires(Component const& component, std::vector<SolverOutputType> const& solver_output,
+             requires(Component const& component, MathOutput<std::vector<SolverOutputType>> const& math_output,
                       Idx2DBranch3 const& math_id) {
                  {
-                     output_result<Component>(component, solver_output, math_id)
+                     output_result<Component>(component, math_output, math_id)
                  } -> detail::assignable_to<std::ranges::range_reference_t<ComponentOutput>>;
              }
 constexpr void output_result(MainModelState<ComponentContainer> const& state,
                              MathOutput<std::vector<SolverOutputType>> const& math_output, ComponentOutput output) {
     detail::produce_output<Component, Idx2DBranch3>(
         state, output, [&math_output](Component const& component, Idx2DBranch3 const& math_id) {
-            return output_result<Component>(component, math_output.solver_output, math_id);
+            return output_result<Component>(component, math_output, math_id);
         });
 }
 template <std::derived_from<Base> Component, class ComponentContainer, typename SolverOutputType,
@@ -605,22 +614,6 @@ constexpr void output_result(MainModelState<ComponentContainer> const& state,
             return output_result<Component, ComponentContainer>(component, state, math_output, topo_id);
         });
 }
-template <std::derived_from<Base> Component, class ComponentContainer, typename SolverOutputType,
-          non_owning_view_c ComponentOutput>
-    requires model_component_state_c<MainModelState, ComponentContainer, Component> &&
-             requires(Component const& component, MainModelState<ComponentContainer> const& state,
-                      MathOutput<SolverOutputType> const& math_output, Idx const obj_seq) {
-                 {
-                     output_result<Component>(component, state, math_output, obj_seq)
-                 } -> detail::assignable_to<std::ranges::range_reference_t<ComponentOutput>>;
-             }
-constexpr void output_result(MainModelState<ComponentContainer> const& state,
-                             MathOutput<SolverOutputType> const& math_output, ComponentOutput output) {
-    detail::produce_output<Component, Idx>(
-        state, output, [&state, &math_output](Component const& component, Idx const obj_seq) {
-            return output_result<Component, ComponentContainer>(component, state, math_output, obj_seq);
-        });
-}
 // vector overload
 template <std::derived_from<Base> Component, class ComponentContainer, typename SolverOutputType, class T>
     requires model_component_state_c<MainModelState, ComponentContainer, Component>
diff --git a/tests/cpp_unit_tests/main_core/test_main_core_output.cpp b/tests/cpp_unit_tests/main_core/test_main_core_output.cpp
index bccfb3a94b..b0c5ca29d6 100644
--- a/tests/cpp_unit_tests/main_core/test_main_core_output.cpp
+++ b/tests/cpp_unit_tests/main_core/test_main_core_output.cpp
@@ -220,6 +220,13 @@ template <symmetry_tag sym> void check_null_current_sensor_output(CurrentSensorO
     check_zero<sym>(output.i_angle_residual);
 }
 
+template <typename SolverOutputType> auto make_math_output(std::vector<SolverOutputType> solver_output) {
+    return MathOutput{.solver_output = std::move(solver_output), .optimizer_output = {}, .supernode_output = {}};
+}
+template <typename SolverOutputType> auto make_math_output(SolverOutputType solver_output) {
+    return make_math_output(std::vector{std::move(solver_output)});
+}
+
 } // namespace
 
 TEST_CASE("Test main core output") {
@@ -384,11 +391,12 @@ TEST_CASE_TEMPLATE("Test main core power sensor output energized state", sym, sy
                                      .q_measured = 2.0}};
 
         // An excluded measurement remains energized while its measured object is active.
-        CHECK(output_result(sensor, state, solver_output, object_sequence).energized == 1);
+        CHECK(output_result(sensor, state, make_math_output(std::move(solver_output)), object_sequence).energized == 1);
 
         appliance.set_status(0);
 
-        check_null_power_sensor_output(output_result(sensor, state, solver_output, object_sequence));
+        check_null_power_sensor_output(
+            output_result(sensor, state, make_math_output(std::move(solver_output)), object_sequence));
     };
 
     check_appliance(MeasuredTerminalType::source, 10, 0, get_component<Source>(state.components, ID{10}));
@@ -409,7 +417,7 @@ TEST_CASE_TEMPLATE("Test main core power sensor output energized state", sym, sy
                                      .p_measured = 1.0,
                                      .q_measured = 2.0}};
 
-        auto const output = output_result(sensor, state, solver_output, 0);
+        auto const output = output_result(sensor, state, make_math_output(std::move(solver_output)), 0);
         if (expected_energized) {
             CHECK(output.energized == 1);
         } else {
@@ -447,7 +455,7 @@ TEST_CASE_TEMPLATE("Test main core power sensor output energized state", sym, sy
                                      .p_measured = 1.0,
                                      .q_measured = 2.0}};
 
-        check_null_power_sensor_output(output_result(sensor, state, solver_output, 0));
+        check_null_power_sensor_output(output_result(sensor, state, make_math_output(std::move(solver_output)), 0));
     }
 }
 
@@ -489,7 +497,7 @@ TEST_CASE_TEMPLATE("Test main core power sensor output with reduced component co
                                  .p_measured = 1.0,
                                  .q_measured = 2.0}};
 
-    CHECK(output_result(sensor, state, solver_output, 0).energized == 1);
+    CHECK(output_result(sensor, state, make_math_output(std::move(solver_output)), 0).energized == 1);
 }
 
 TEST_CASE_TEMPLATE("Test main core current sensor output energized state", sym, symmetric_t, asymmetric_t) {
@@ -511,7 +519,7 @@ TEST_CASE_TEMPLATE("Test main core current sensor output energized state", sym,
                                        .i_angle_measured = 0.0},
                                       u_rated};
 
-        auto const output = output_result(sensor, state, solver_output, 0);
+        auto const output = output_result(sensor, state, make_math_output(std::move(solver_output)), 0);
         if (terminal_type == MeasuredTerminalType::branch_from) {
             check_null_current_sensor_output(output);
         } else {
@@ -529,7 +537,8 @@ TEST_CASE_TEMPLATE("Test main core current sensor output energized state", sym,
                                            .i_measured = 1.0,
                                            .i_angle_measured = 0.0},
                                           u_rated};
-    check_null_current_sensor_output(output_result(branch3_sensor, state, solver_output, 0));
+    check_null_current_sensor_output(
+        output_result(branch3_sensor, state, make_math_output(std::move(solver_output)), 0));
 
     auto disconnected_coupling = std::make_shared<TopologicalComponentToMathCoupling>(*state.topo_comp_coup);
     disconnected_coupling->branch = {{.group = disconnected, .pos = disconnected}};
@@ -555,7 +564,7 @@ TEST_CASE_TEMPLATE("Test main core current sensor output energized state", sym,
                                        .i_angle_measured = 0.0},
                                       u_rated};
 
-        check_null_current_sensor_output(output_result(sensor, state, solver_output, 0));
+        check_null_current_sensor_output(output_result(sensor, state, make_math_output(std::move(solver_output)), 0));
     }
 }
 } // namespace power_grid_model::main_core