diff --git a/include/scalfmm/matrix_kernels/gaussian.hpp b/include/scalfmm/matrix_kernels/gaussian.hpp
index f1e548f9dbfb39e2338f82c8bfcc048162904071..127501e7dda557193c76cb80deae609ef675ca66 100644
--- a/include/scalfmm/matrix_kernels/gaussian.hpp
+++ b/include/scalfmm/matrix_kernels/gaussian.hpp
@@ -40,7 +40,8 @@ namespace scalfmm::matrix_kernels
         // Optional constant
         static constexpr bool is_smooth = true;   // Specify the kernel is smooth and can be evaluated at x == y.
 
-        value_type m_coeff{value_type(1.)};   // parameter/coefficient of the kernel.
+        value_type m_coeff{value_type(1.)};     // parameter/coefficient of the kernel.
+        value_type m_epsilon{value_type(0.)};   // epsilon value for the evaluation at point x (self-interaction).
 
         /**
          * @brief Set the value of the coefficient of the kernel.
@@ -49,12 +50,23 @@ namespace scalfmm::matrix_kernels
          */
         void set_coeff(value_type in_coeff) { m_coeff = in_coeff; }
 
+        /**
+         * @brief Set the value of epsilon (to compute self-interaction).
+         *
+         * @param in_epsilon value of epsilon to be set.
+         */
+        void set_epsilon(value_type in_epsilon) { m_epsilon = in_epsilon; }
+
         /**
 	* @brief Returns the name of the kernel.
 	*
 	* @return A string representing the kernel's name.
 	*/
-        const std::string name() const { return std::string("gaussian ") + "  coeff = " + std::to_string(m_coeff); }
+        const std::string name() const
+        {
+            return std::string("gaussian") + ", coeff = " + std::to_string(m_coeff) +
+                   ", m_epsilon = " + std::to_string(m_epsilon);
+        }
 
         /**
 	* @brief Returns the mutual coefficient of size \f$kn\f$.
@@ -72,6 +84,43 @@ namespace scalfmm::matrix_kernels
             return vector_type<ValueType>{ValueType(1.)};
         }
 
+        /**
+	 * @brief Overload of the `()` operator to evaluate the kernel at point \f$x\f$ (self-interaction).
+	*
+	* @tparam ValueType Type of the elements in the point.
+	* @tparam Dimension The dimension of the point.
+	*
+	* @param x Multidimensional point (e.g., a 2D point).
+	*
+	* @return The matrix \f$K(x, x)\f$ in a vector (row-major storage).
+	*/
+        template<typename ValueType, std::size_t Dimension>
+        [[nodiscard]] inline auto operator()(container::point<ValueType, Dimension> const& x) const noexcept
+          -> matrix_type<std::decay_t<ValueType>>
+        {
+            return evaluate(x);
+        }
+
+        /**
+	 * @brief Evaluates the kernel at points \f$x\f$ (self-interaction).
+	*
+	* @tparam ValueType Type of the elements in the point.
+	* @tparam Dimension The dimension of the point.
+	*
+	* @param x Multidimensional point (e.g., a 2D point).
+	*
+	* @return The matrix \f$K(x, x)\f$ in a vector (row-major storage).
+	*/
+        template<typename ValueType, std::size_t Dimension>
+        [[nodiscard]] inline auto evaluate(container::point<ValueType, Dimension> const& x) const noexcept
+          -> matrix_type<std::decay_t<ValueType>>
+        {
+            using decayed_type = std::decay_t<ValueType>;
+            matrix_type<decayed_type> tmp;
+            std::fill(std::begin(tmp), std::end(tmp), decayed_type(1. + m_epsilon));
+            return tmp;
+        }
+
         /**
 	* @brief Overload of the `()` operator to evaluate the kernel at points \f$x\f$ and \f$y\f$.
 	*