diff --git a/framework/include/physics/PhysicsBase.h b/framework/include/physics/PhysicsBase.h
index 4bee0d914884..93213f0bd617 100644
--- a/framework/include/physics/PhysicsBase.h
+++ b/framework/include/physics/PhysicsBase.h
@@ -94,8 +94,8 @@ class PhysicsBase : public Action, public InputParametersChecksUtils<PhysicsBase
   /// restriction of the Physics. More complex behavior should be implemented by overriding
   virtual void addComponent(const ActionComponent & component);
 
-  /// Return the list of nonlinear variables in this physics
-  const std::vector<VariableName> & nonlinearVariableNames() const { return _nl_var_names; };
+  /// Return the list of solver (nonlinear + linear) variables in this physics
+  const std::vector<VariableName> & solverVariableNames() const { return _solver_var_names; };
   /// Return the list of aux variables in this physics
   const std::vector<VariableName> & auxVariableNames() const { return _aux_var_names; };
 
@@ -134,14 +134,23 @@ class PhysicsBase : public Action, public InputParametersChecksUtils<PhysicsBase
   /// Use prefix() to disambiguate names
   std::string prefix() const { return name() + "_"; }
 
-  /// Keep track of the name of a nonlinear variable defined in the Physics
-  void saveSolverVariableName(const VariableName & var_name) { _nl_var_names.push_back(var_name); }
+  /// Keep track of the name of the solver variable defined in the Physics
+  void saveSolverVariableName(const VariableName & var_name)
+  {
+    _solver_var_names.push_back(var_name);
+  }
   /// Keep track of the name of an aux variable defined in the Physics
   void saveAuxVariableName(const VariableName & var_name) { _aux_var_names.push_back(var_name); }
 
   /// Check whether a variable already exists
   bool variableExists(const VariableName & var_name, bool error_if_aux) const;
 
+  /// Get the solver system for this variable index. The index should be the index of the variable in solver
+  /// var_names (currently _solver_var_names) vector
+  const SolverSystemName & getSolverSystem(unsigned int variable_index) const;
+  /// Get the solver system for this variable name
+  const SolverSystemName & getSolverSystem(const VariableName & variable_name) const;
+
   /// Add a new required task for all physics deriving from this class
   /// NOTE: This does not register the task, you still need to call registerMooseAction
   void addRequiredPhysicsTask(const std::string & task) { _required_tasks.insert(task); }
@@ -158,8 +167,11 @@ class PhysicsBase : public Action, public InputParametersChecksUtils<PhysicsBase
    */
   bool allMeshBlocks(const std::vector<SubdomainName> & blocks) const;
 
-  /// System number for the system owning the variables
-  const unsigned int _sys_number;
+  /// System names for the system(s) owning the solver variables
+  std::vector<SolverSystemName> _system_names;
+
+  /// System numbers for the system(s) owning the solver variables
+  std::vector<unsigned int> _system_numbers;
 
   /// Whether to output additional information
   const bool _verbose;
@@ -224,8 +236,8 @@ class PhysicsBase : public Action, public InputParametersChecksUtils<PhysicsBase
   /// some physics directly to steady state
   MooseEnum _is_transient;
 
-  /// Vector of the nonlinear variables in the Physics
-  std::vector<VariableName> _nl_var_names;
+  /// Vector of the solver variables (nonlinear and linear) in the Physics
+  std::vector<VariableName> _solver_var_names;
   /// Vector of the aux variables in the Physics
   std::vector<VariableName> _aux_var_names;
 
diff --git a/framework/src/physics/PhysicsBase.C b/framework/src/physics/PhysicsBase.C
index bfd222af7f6d..e7fbad5f3e05 100644
--- a/framework/src/physics/PhysicsBase.C
+++ b/framework/src/physics/PhysicsBase.C
@@ -29,6 +29,12 @@ PhysicsBase::validParams()
   params.addParam<MooseEnum>(
       "transient", transient_options, "Whether the physics is to be solved as a transient");
 
+  params.addParam<std::vector<SolverSystemName>>(
+      "system_name",
+      {"nl0"},
+      "Name of the solver system(s) for the variables. If a single name is specified, "
+      "that system is used for all solver variables.");
+
   MooseEnum pc_options("default none", "none");
   params.addParam<MooseEnum>(
       "preconditioning", pc_options, "Which preconditioning to use for this Physics");
@@ -53,7 +59,7 @@ PhysicsBase::validParams()
 PhysicsBase::PhysicsBase(const InputParameters & parameters)
   : Action(parameters),
     InputParametersChecksUtils<PhysicsBase>(this),
-    _sys_number(0),
+    _system_names(getParam<std::vector<SolverSystemName>>("system_name")),
     _verbose(getParam<bool>("verbose")),
     _preconditioning(getParam<MooseEnum>("preconditioning")),
     _blocks(getParam<std::vector<SubdomainName>>("block")),
@@ -143,7 +149,7 @@ PhysicsBase::act()
   // Exodus restart capabilities
   if (_current_task == "copy_vars_physics")
   {
-    copyVariablesFromMesh(nonlinearVariableNames(), true);
+    copyVariablesFromMesh(solverVariableNames(), true);
     if (_aux_var_names.size() > 0)
       copyVariablesFromMesh(auxVariableNames(), false);
   }
@@ -244,6 +250,19 @@ PhysicsBase::initializePhysics()
   else
     _dim = _mesh->dimension();
 
+  // Check that the systems exist in the Problem
+  // TODO: try to add the systems to the problem from here
+  const auto & problem_systems = getProblem().getNonlinearSystemNames();
+  for (const auto & sys_name : _system_names)
+    if (std::find(problem_systems.begin(), problem_systems.end(), sys_name) ==
+            problem_systems.end() &&
+        solverVariableNames().size())
+      mooseError("System '", sys_name, "' is not used in the problem");
+
+  // Cache system number as some routines prefer those
+  for (const auto & sys_name : _system_names)
+    _system_numbers.push_back(getProblem().solverSysNum(sys_name));
+
   // Forward physics verbosity to problem to output the setup
   if (_verbose)
     getProblem().setVerboseProblem(_verbose);
@@ -257,6 +276,14 @@ PhysicsBase::checkIntegrityEarly() const
 {
   if (_is_transient == "true" && !getProblem().isTransient())
     paramError("transient", "We cannot solve a physics as transient in a steady problem");
+
+  // Check that there is a system for each variable
+  if (_system_names.size() != 1 && _system_names.size() != _solver_var_names.size())
+    paramError("system_name",
+               "There should be one system name per solver variable (potentially repeated), or a "
+               "single system name for all variables. Current you have '" +
+                   std::to_string(_system_names.size()) + "' systems specified for '" +
+                   std::to_string(_solver_var_names.size()) + "' solver variables.");
 }
 
 void
@@ -265,14 +292,19 @@ PhysicsBase::copyVariablesFromMesh(const std::vector<VariableName> & variables_t
 {
   if (getParam<bool>("initialize_variables_from_mesh_file"))
   {
-    SystemBase & system = are_nonlinear ? getProblem().getNonlinearSystemBase(_sys_number)
-                                        : getProblem().systemBaseAuxiliary();
     mooseInfoRepeated("Adding Exodus restart for " + std::to_string(variables_to_copy.size()) +
                       " variables: " + Moose::stringify(variables_to_copy));
     // TODO Check that the variable types and orders are actually supported for exodus restart
-    for (const auto & var_name : variables_to_copy)
+    for (const auto i : index_range(variables_to_copy))
+    {
+      SystemBase & system = are_nonlinear
+                                ? getProblem().getNonlinearSystemBase(
+                                      _system_numbers.size() == 1 ? _system_numbers[0] : i)
+                                : getProblem().systemBaseAuxiliary();
+      const auto & var_name = variables_to_copy[i];
       system.addVariableToCopy(
           var_name, var_name, getParam<std::string>("initial_from_file_timestep"));
+    }
   }
 }
 
@@ -291,6 +323,31 @@ PhysicsBase::variableExists(const VariableName & var_name, bool error_if_aux) co
     return false;
 }
 
+const SolverSystemName &
+PhysicsBase::getSolverSystem(unsigned int variable_index) const
+{
+  mooseAssert(!_system_names.empty(), "We should have a solver system name");
+  if (_system_names.size() == 1)
+    return _system_names[0];
+  else
+    return _system_names[variable_index];
+}
+
+const SolverSystemName &
+PhysicsBase::getSolverSystem(const VariableName & var_name) const
+{
+  mooseAssert(!_system_names.empty(), "We should have a solver system name");
+  // No need to look if only one system
+  if (_system_names.size() == 1)
+    return _system_names[0];
+
+  // Find the variable in the list of variables
+  for (const auto variable_index : index_range(_solver_var_names))
+    if (var_name == _solver_var_names[variable_index])
+      return _system_names[variable_index];
+  mooseError("Variable ", var_name, " was not found in list of variables");
+}
+
 void
 PhysicsBase::checkRequiredTasks() const
 {
diff --git a/framework/src/problems/FEProblemBase.C b/framework/src/problems/FEProblemBase.C
index d7031b1783d8..d5e2b67eec10 100644
--- a/framework/src/problems/FEProblemBase.C
+++ b/framework/src/problems/FEProblemBase.C
@@ -6053,7 +6053,7 @@ FEProblemBase::init()
     _aux->dofMap().attach_extra_send_list_function(&extraSendList, _aux.get());
 
     if (!_skip_nl_system_check && _solve && n_vars == 0)
-      mooseError("No variables specified in the FEProblemBase '", name(), "'.");
+      mooseError("No variables specified in nonlinear system '", nl->name(), "'.");
   }
 
   ghostGhostedBoundaries(); // We do this again right here in case new boundaries have been added