diff --git a/src/FEM/FEMQuadratureData.h b/src/FEM/FEMQuadratureData.h
new file mode 100644
index 000000000..1824c572a
--- /dev/null
+++ b/src/FEM/FEMQuadratureData.h
@@ -0,0 +1,25 @@
+#ifndef IPPL_FEM_QUADRATURE_DATA_H
+#define IPPL_FEM_QUADRATURE_DATA_H
+
+#include "Types/Vector.h"
+
+namespace ippl {
+
+/**
+ * @brief Per-quadrature-node basis data passed to evaluateAx evaluator functors.
+ *
+ * @tparam TVal Type of basis values at each local DOF (scalar T for Lagrange,
+ *              Vector<T, Dim> for Nedelec).
+ * @tparam TDeriv Type of the spatial derivative data at each local DOF
+ *                (gradient for Lagrange, curl for Nedelec).
+ * @tparam numElementDOFs Number of local DOFs per element.
+ */
+template <typename TVal, typename TDeriv, unsigned numElementDOFs>
+struct QuadratureData {
+    const Vector<TVal, numElementDOFs>& val_q;
+    const Vector<TDeriv, numElementDOFs>& deriv_q;
+};
+
+}  // namespace ippl
+
+#endif
diff --git a/src/FEM/LagrangeSpace.h b/src/FEM/LagrangeSpace.h
index dcbfdd559..b744cd49b 100644
--- a/src/FEM/LagrangeSpace.h
+++ b/src/FEM/LagrangeSpace.h
@@ -7,6 +7,7 @@
 
 #include <cmath>
 
+#include "FEM/FEMQuadratureData.h"
 #include "FEM/FiniteElementSpace.h"
 
 constexpr unsigned getLagrangeNumElementDOFs(unsigned Dim, unsigned Order) {
diff --git a/src/FEM/LagrangeSpace.hpp b/src/FEM/LagrangeSpace.hpp
index b1743cd4f..f13385014 100644
--- a/src/FEM/LagrangeSpace.hpp
+++ b/src/FEM/LagrangeSpace.hpp
@@ -394,11 +394,12 @@ namespace ippl {
         const Vector<point_t, QuadratureType::numElementNodes> q =
             this->quadrature_m.getIntegrationNodesForRefElement();
 
-        // TODO move outside of evaluateAx (I think it is possible for other problems as well)
-        // Gradients of the basis functions for the DOF at the quadrature nodes
+        // Basis function values and gradients at the quadrature nodes
+        Vector<Vector<T, numElementDOFs>, QuadratureType::numElementNodes> b_q;
         Vector<Vector<point_t, numElementDOFs>, QuadratureType::numElementNodes> grad_b_q;
         for (size_t k = 0; k < QuadratureType::numElementNodes; ++k) {
             for (size_t i = 0; i < numElementDOFs; ++i) {
+                b_q[k][i]      = this->evaluateRefElementShapeFunction(i, q[k]);
                 grad_b_q[k][i] = this->evaluateRefElementShapeFunctionGradient(i, q[k]);
             }
         }
@@ -412,7 +413,8 @@ namespace ippl {
             for (size_t j = 0; j < numElementDOFs; ++j) {
                 A_K[i][j] = 0.0;
                 for (size_t k = 0; k < QuadratureType::numElementNodes; ++k) {
-                    A_K[i][j] += w[k] * evalFunction(i, j, grad_b_q[k]);
+                    A_K[i][j] += w[k] * evalFunction(
+                        i, j, QuadratureData<T, point_t, numElementDOFs>{b_q[k], grad_b_q[k]});
                 }
             }
         }
@@ -544,11 +546,12 @@ namespace ippl {
         const Vector<point_t, QuadratureType::numElementNodes> q =
             this->quadrature_m.getIntegrationNodesForRefElement();
 
-        // TODO move outside of evaluateAx (I think it is possible for other problems as well)
-        // Gradients of the basis functions for the DOF at the quadrature nodes
+        // Basis function values and gradients at the quadrature nodes
+        Vector<Vector<T, numElementDOFs>, QuadratureType::numElementNodes> b_q;
         Vector<Vector<point_t, numElementDOFs>, QuadratureType::numElementNodes> grad_b_q;
         for (size_t k = 0; k < QuadratureType::numElementNodes; ++k) {
             for (size_t i = 0; i < numElementDOFs; ++i) {
+                b_q[k][i]      = this->evaluateRefElementShapeFunction(i, q[k]);
                 grad_b_q[k][i] = this->evaluateRefElementShapeFunctionGradient(i, q[k]);
             }
         }
@@ -562,7 +565,8 @@ namespace ippl {
             for (size_t j = 0; j < numElementDOFs; ++j) {
                 A_K[i][j] = 0.0;
                 for (size_t k = 0; k < QuadratureType::numElementNodes; ++k) {
-                    A_K[i][j] += w[k] * evalFunction(i, j, grad_b_q[k]);
+                    A_K[i][j] += w[k] * evalFunction(
+                        i, j, QuadratureData<T, point_t, numElementDOFs>{b_q[k], grad_b_q[k]});
                 }
             }
         }
@@ -686,11 +690,12 @@ namespace ippl {
         const Vector<point_t, QuadratureType::numElementNodes> q =
             this->quadrature_m.getIntegrationNodesForRefElement();
 
-        // TODO move outside of evaluateAx (I think it is possible for other problems as well)
-        // Gradients of the basis functions for the DOF at the quadrature nodes
+        // Basis function values and gradients at the quadrature nodes
+        Vector<Vector<T, numElementDOFs>, QuadratureType::numElementNodes> b_q;
         Vector<Vector<point_t, numElementDOFs>, QuadratureType::numElementNodes> grad_b_q;
         for (size_t k = 0; k < QuadratureType::numElementNodes; ++k) {
             for (size_t i = 0; i < numElementDOFs; ++i) {
+                b_q[k][i]      = this->evaluateRefElementShapeFunction(i, q[k]);
                 grad_b_q[k][i] = this->evaluateRefElementShapeFunctionGradient(i, q[k]);
             }
         }
@@ -704,7 +709,8 @@ namespace ippl {
             for (size_t j = 0; j < numElementDOFs; ++j) {
                 A_K[i][j] = 0.0;
                 for (size_t k = 0; k < QuadratureType::numElementNodes; ++k) {
-                    A_K[i][j] += w[k] * evalFunction(i, j, grad_b_q[k]);
+                    A_K[i][j] += w[k] * evalFunction(
+                        i, j, QuadratureData<T, point_t, numElementDOFs>{b_q[k], grad_b_q[k]});
                 }
             }
         }
@@ -829,10 +835,12 @@ namespace ippl {
             this->quadrature_m.getIntegrationNodesForRefElement();
 
         // TODO move outside of evaluateAx (I think it is possible for other problems as well)
-        // Gradients of the basis functions for the DOF at the quadrature nodes
+        // Basis function values and gradients at the quadrature nodes
+        Vector<Vector<T, numElementDOFs>, QuadratureType::numElementNodes> b_q;
         Vector<Vector<point_t, numElementDOFs>, QuadratureType::numElementNodes> grad_b_q;
         for (size_t k = 0; k < QuadratureType::numElementNodes; ++k) {
             for (size_t i = 0; i < numElementDOFs; ++i) {
+                b_q[k][i]      = this->evaluateRefElementShapeFunction(i, q[k]);
                 grad_b_q[k][i] = this->evaluateRefElementShapeFunctionGradient(i, q[k]);
             }
         }
@@ -846,7 +854,8 @@ namespace ippl {
             for (size_t j = 0; j < numElementDOFs; ++j) {
                 A_K[i][j] = 0.0;
                 for (size_t k = 0; k < QuadratureType::numElementNodes; ++k) {
-                    A_K[i][j] += w[k] * evalFunction(i, j, grad_b_q[k]);
+                    A_K[i][j] += w[k] * evalFunction(
+                        i, j, QuadratureData<T, point_t, numElementDOFs>{b_q[k], grad_b_q[k]});
                 }
             }
         }
@@ -967,11 +976,12 @@ namespace ippl {
         const Vector<point_t, QuadratureType::numElementNodes> q =
             this->quadrature_m.getIntegrationNodesForRefElement();
 
-        // TODO move outside of evaluateAx (I think it is possible for other problems as well)
-        // Gradients of the basis functions for the DOF at the quadrature nodes
+        // Basis function values and gradients at the quadrature nodes
+        Vector<Vector<T, numElementDOFs>, QuadratureType::numElementNodes> b_q;
         Vector<Vector<point_t, numElementDOFs>, QuadratureType::numElementNodes> grad_b_q;
         for (size_t k = 0; k < QuadratureType::numElementNodes; ++k) {
             for (size_t i = 0; i < numElementDOFs; ++i) {
+                b_q[k][i]      = this->evaluateRefElementShapeFunction(i, q[k]);
                 grad_b_q[k][i] = this->evaluateRefElementShapeFunctionGradient(i, q[k]);
             }
         }
@@ -984,7 +994,8 @@ namespace ippl {
         for (size_t i = 0; i < numElementDOFs; ++i) {
             A_K_diag[i] = 0.0;
             for (size_t k = 0; k < QuadratureType::numElementNodes; ++k) {
-                A_K_diag[i] += w[k] * evalFunction(i, i, grad_b_q[k]);
+                A_K_diag[i] += w[k] * evalFunction(
+                    i, i, QuadratureData<T, point_t, numElementDOFs>{b_q[k], grad_b_q[k]});
             }
         }
 
@@ -1096,11 +1107,12 @@ namespace ippl {
         const Vector<point_t, QuadratureType::numElementNodes> q =
             this->quadrature_m.getIntegrationNodesForRefElement();
 
-        // TODO move outside of evaluateAx (I think it is possible for other problems as well)
-        // Gradients of the basis functions for the DOF at the quadrature nodes
+        // Basis function values and gradients at the quadrature nodes
+        Vector<Vector<T, numElementDOFs>, QuadratureType::numElementNodes> b_q;
         Vector<Vector<point_t, numElementDOFs>, QuadratureType::numElementNodes> grad_b_q;
         for (size_t k = 0; k < QuadratureType::numElementNodes; ++k) {
             for (size_t i = 0; i < numElementDOFs; ++i) {
+                b_q[k][i]      = this->evaluateRefElementShapeFunction(i, q[k]);
                 grad_b_q[k][i] = this->evaluateRefElementShapeFunctionGradient(i, q[k]);
             }
         }
@@ -1113,7 +1125,8 @@ namespace ippl {
         for (size_t i = 0; i < numElementDOFs; ++i) {
             A_K_diag[i] = 0.0;
             for (size_t k = 0; k < QuadratureType::numElementNodes; ++k) {
-                A_K_diag[i] += w[k] * evalFunction(i, i, grad_b_q[k]);
+                A_K_diag[i] += w[k] * evalFunction(
+                    i, i, QuadratureData<T, point_t, numElementDOFs>{b_q[k], grad_b_q[k]});
             }
         }
 
@@ -1216,11 +1229,12 @@ namespace ippl {
         const Vector<point_t, QuadratureType::numElementNodes> q =
             this->quadrature_m.getIntegrationNodesForRefElement();
 
-        // TODO move outside of evaluateAx (I think it is possible for other problems as well)
-        // Gradients of the basis functions for the DOF at the quadrature nodes
+        // Basis function values and gradients at the quadrature nodes
+        Vector<Vector<T, numElementDOFs>, QuadratureType::numElementNodes> b_q;
         Vector<Vector<point_t, numElementDOFs>, QuadratureType::numElementNodes> grad_b_q;
         for (size_t k = 0; k < QuadratureType::numElementNodes; ++k) {
             for (size_t i = 0; i < numElementDOFs; ++i) {
+                b_q[k][i]      = this->evaluateRefElementShapeFunction(i, q[k]);
                 grad_b_q[k][i] = this->evaluateRefElementShapeFunctionGradient(i, q[k]);
             }
         }
@@ -1234,7 +1248,8 @@ namespace ippl {
             for (size_t j = 0; j < numElementDOFs; ++j) {
                 A_K[i][j] = 0.0;
                 for (size_t k = 0; k < QuadratureType::numElementNodes; ++k) {
-                    A_K[i][j] += w[k] * evalFunction(i, j, grad_b_q[k]);
+                    A_K[i][j] += w[k] * evalFunction(
+                        i, j, QuadratureData<T, point_t, numElementDOFs>{b_q[k], grad_b_q[k]});
                 }
             }
         }
diff --git a/src/FEM/NedelecSpace.h b/src/FEM/NedelecSpace.h
index c023cbf34..f99745c3d 100644
--- a/src/FEM/NedelecSpace.h
+++ b/src/FEM/NedelecSpace.h
@@ -7,6 +7,7 @@
 
 #include <cmath>
 
+#include "FEM/FEMQuadratureData.h"
 #include "FEM/FEMVector.h"
 #include "FEM/FiniteElementSpace.h"
 
diff --git a/src/FEM/NedelecSpace.hpp b/src/FEM/NedelecSpace.hpp
index bb3888356..5fefc0f1f 100644
--- a/src/FEM/NedelecSpace.hpp
+++ b/src/FEM/NedelecSpace.hpp
@@ -425,7 +425,10 @@ namespace ippl {
             for (size_t j = 0; j < numElementDOFs; ++j) {
                 A[i][j] = 0.0;
                 for (size_t k = 0; k < QuadratureType::numElementNodes; ++k) {
-                    A[i][j] += w[k] * evalFunction(i, j, curl_b_q[k], val_b_q[k]);
+                    A[i][j] += w[k] * evalFunction(
+                        i, j,
+                        QuadratureData<Vector<T, Dim>, Vector<T, Dim>, numElementDOFs>{val_b_q[k],
+                                                                                      curl_b_q[k]});
                 }
             }
         }
diff --git a/src/MaxwellSolvers/FEMMaxwellDiffusionSolver.h b/src/MaxwellSolvers/FEMMaxwellDiffusionSolver.h
index 0319ce8c4..82c4b5f9b 100644
--- a/src/MaxwellSolvers/FEMMaxwellDiffusionSolver.h
+++ b/src/MaxwellSolvers/FEMMaxwellDiffusionSolver.h
@@ -9,6 +9,7 @@
 #include <iomanip>
 #include <iostream>
 
+#include "FEM/FEMQuadratureData.h"
 #include "LinearSolvers/PCG.h"
 #include "Maxwell.h"
 
@@ -62,17 +63,15 @@ namespace ippl {
          *
          * @param i The first DOF index.
          * @param j The second DOF index.
-         * @param curl_b_q_k The curl of the DOFs.
-         * @param val_b_q_k The values of the DOFs.
+         * @param qd Per-quadrature basis values and curls.
          *
          * @returns (curl(b_i)*curl(b_j) + b_i*b_j)*absDetDPhi
          */
         KOKKOS_FUNCTION auto operator()(
             size_t i, size_t j,
-            const ippl::Vector<ippl::Vector<T, Dim>, numElementDOFs>& curl_b_q_k,
-            const ippl::Vector<ippl::Vector<T, Dim>, numElementDOFs>& val_b_q_k) const {
-            T curlTerm = dot(DPhiInvT * curl_b_q_k[j], DPhiInvT * curl_b_q_k[i]).apply();
-            T massTerm = dot(val_b_q_k[j], val_b_q_k[i]).apply();
+            const QuadratureData<Vector<T, Dim>, Vector<T, Dim>, numElementDOFs>& qd) const {
+            T curlTerm = dot(DPhiInvT * qd.deriv_q[j], DPhiInvT * qd.deriv_q[i]).apply();
+            T massTerm = dot(qd.val_q[j], qd.val_q[i]).apply();
             return (curlTerm + massTerm) * absDetDPhi;
         }
     };
diff --git a/src/PoissonSolvers/EvalFunctor.h b/src/PoissonSolvers/EvalFunctor.h
index f757f1cc6..ce9918ba7 100644
--- a/src/PoissonSolvers/EvalFunctor.h
+++ b/src/PoissonSolvers/EvalFunctor.h
@@ -9,6 +9,8 @@
 #ifndef IPPL_EVALFUNCTOR_H
 #define IPPL_EVALFUNCTOR_H
 
+#include "FEM/FEMQuadratureData.h"
+
 namespace ippl {
     template <typename Tlhs, unsigned Dim, unsigned numElemDOFs>
     struct EvalFunctor {
@@ -21,8 +23,8 @@ namespace ippl {
 
         KOKKOS_FUNCTION auto operator()(
             const size_t& i, const size_t& j,
-            const Vector<Vector<Tlhs, Dim>, numElemDOFs>& grad_b_q_k) const {
-            return dot((DPhiInvT * grad_b_q_k[j]), (DPhiInvT * grad_b_q_k[i])).apply() * absDetDPhi;
+            const QuadratureData<Tlhs, Vector<Tlhs, Dim>, numElemDOFs>& qd) const {
+            return dot((DPhiInvT * qd.deriv_q[j]), (DPhiInvT * qd.deriv_q[i])).apply() * absDetDPhi;
         }
     };
 }
diff --git a/unit_tests/FEM/CMakeLists.txt b/unit_tests/FEM/CMakeLists.txt
index b15a1273f..ab1b5456c 100644
--- a/unit_tests/FEM/CMakeLists.txt
+++ b/unit_tests/FEM/CMakeLists.txt
@@ -7,6 +7,7 @@ add_ippl_test(EdgeElement)
 add_ippl_test(QuadrilateralElement)
 add_ippl_test(HexahedralElement)
 add_ippl_test(FiniteElementSpace)
+add_ippl_test(FEMQuadratureData)
 add_ippl_test(LagrangeSpace)
 add_ippl_test(Quadrature)
 add_ippl_test(GaussJacobiQuadrature)
diff --git a/unit_tests/FEM/FEMQuadratureData.cpp b/unit_tests/FEM/FEMQuadratureData.cpp
new file mode 100644
index 000000000..e9ae7835a
--- /dev/null
+++ b/unit_tests/FEM/FEMQuadratureData.cpp
@@ -0,0 +1,44 @@
+#include "Ippl.h"
+
+#include "FEM/FEMQuadratureData.h"
+
+#include "gtest/gtest.h"
+
+TEST(FEMQuadratureDataTest, DataAccess) {
+    using T = double;
+    constexpr unsigned Dim       = 2;
+    constexpr unsigned numDOFs   = 4;
+    using point_t                = ippl::Vector<T, Dim>;
+
+    ippl::Vector<T, numDOFs> val;
+    val[0] = 0.1;
+    val[1] = 0.2;
+    val[2] = 0.3;
+    val[3] = 0.4;
+
+    ippl::Vector<point_t, numDOFs> deriv;
+    deriv[0] = point_t(1.0, 0.0);
+    deriv[1] = point_t(0.0, 1.0);
+    deriv[2] = point_t(-1.0, 2.0);
+    deriv[3] = point_t(0.5, -0.5);
+
+    const ippl::QuadratureData<T, point_t, numDOFs> qd{val, deriv};
+
+    EXPECT_DOUBLE_EQ(qd.val_q[0], 0.1);
+    EXPECT_DOUBLE_EQ(qd.val_q[2], 0.3);
+    EXPECT_DOUBLE_EQ(qd.deriv_q[1](0), 0.0);
+    EXPECT_DOUBLE_EQ(qd.deriv_q[1](1), 1.0);
+    EXPECT_DOUBLE_EQ(qd.deriv_q[3](0), 0.5);
+    EXPECT_DOUBLE_EQ(qd.deriv_q[3](1), -0.5);
+
+    val[2] = 0.35;
+    deriv[0] = point_t(2.0, 3.0);
+    EXPECT_DOUBLE_EQ(qd.val_q[2], 0.35);
+    EXPECT_DOUBLE_EQ(qd.deriv_q[0](0), 2.0);
+    EXPECT_DOUBLE_EQ(qd.deriv_q[0](1), 3.0);
+}
+
+int main(int argc, char* argv[]) {
+    ::testing::InitGoogleTest(&argc, argv);
+    return RUN_ALL_TESTS();
+}
diff --git a/unit_tests/FEM/LagrangeSpace.cpp b/unit_tests/FEM/LagrangeSpace.cpp
index 79dcfdd96..0adae2553 100644
--- a/unit_tests/FEM/LagrangeSpace.cpp
+++ b/unit_tests/FEM/LagrangeSpace.cpp
@@ -20,8 +20,8 @@ struct EvalFunctor {
         , absDetDPhi(absDetDPhi) {}
 
     KOKKOS_FUNCTION auto operator()(const size_t& i, const size_t& j,
-                    const ippl::Vector<ippl::Vector<Tlhs, Dim>, numElemDOFs>& grad_b_q_k) const {
-        return dot((DPhiInvT * grad_b_q_k[j]), (DPhiInvT * grad_b_q_k[i])).apply() * absDetDPhi;
+                    const ippl::QuadratureData<Tlhs, ippl::Vector<Tlhs, Dim>, numElemDOFs>& qd) const {
+        return dot((DPhiInvT * qd.deriv_q[j]), (DPhiInvT * qd.deriv_q[i])).apply() * absDetDPhi;
     }
 };
 
diff --git a/unit_tests/FEM/NedelecSpace.cpp b/unit_tests/FEM/NedelecSpace.cpp
index b89253954..59db5d789 100644
--- a/unit_tests/FEM/NedelecSpace.cpp
+++ b/unit_tests/FEM/NedelecSpace.cpp
@@ -12,8 +12,8 @@ class NedelecSpaceTest;
 template <typename T, unsigned Dim, unsigned numElementDOFs>
 struct DummyFunctor {
     KOKKOS_FUNCTION const T operator()(size_t i, size_t j,
-        [[maybe_unused]] const ippl::Vector<ippl::Vector<T, Dim>, numElementDOFs>& curl_b_q_k,
-        [[maybe_unused]] const ippl::Vector<ippl::Vector<T, Dim>, numElementDOFs>& val_b_q_k) const {
+        [[maybe_unused]] const ippl::QuadratureData<ippl::Vector<T, Dim>,
+        ippl::Vector<T, Dim>, numElementDOFs>& qd) const {
         return i==j; //val_b_q_k<:i:>.dot(val_b_q_k<:j:>);
     }
 };