diff --git a/kernel_embedding_dictionary/embeddings/embedding.py b/kernel_embedding_dictionary/embeddings/embedding.py
index f55fdb6..9c5ee34 100644
--- a/kernel_embedding_dictionary/embeddings/embedding.py
+++ b/kernel_embedding_dictionary/embeddings/embedding.py
@@ -19,6 +19,9 @@
     matern72_lebesgue_mean_func_1d,
     matern_lebesgue_mean_func_1d,
 )
+from .var_funcs_1d import (
+    expquad_lebesgue_var_func_1d,
+)
 
 
 class KernelEmbedding:
@@ -31,7 +34,8 @@ def __init__(self, kernel: ProductKernel, measure: ProductMeasure):
         self._measure = measure
 
         # kernel and measure must be set first
-        self._mean_func_1d = self._get_1d_funcs()
+        self._mean_func_1d = self._get_1d_mean_funcs()
+        self._var_func_1d = self._get_1d_var_funcs()
 
     def __str__(self) -> str:
         return f"Kernel embedding for\n\n{self._kernel.__str__()}\n\nand\n\n{self._measure.__str__()}"
@@ -52,8 +56,15 @@ def mean(self, x: np.ndarray) -> np.ndarray:
             params_dim = {**self._kernel.get_param_dict_from_dim(dim), **self._measure.get_param_dict_from_dim(dim)}
             kernel_mean *= self._mean_func_1d(x[:, dim], **params_dim)
         return kernel_mean
+    
+    def variance(self) -> float:
+        kernel_var = 1
+        for dim in range(self.ndim):
+            params_dim = {**self._kernel.get_param_dict_from_dim(dim), **self._measure.get_param_dict_from_dim(dim)}
+            kernel_var *= self._var_func_1d(**params_dim)
+        return kernel_var
 
-    def _get_1d_funcs(self) -> Callable:
+    def _get_1d_mean_funcs(self) -> Callable:
 
         mean_func_1d_dict = {
             "expquad-lebesgue": expquad_lebesgue_mean_func_1d,
@@ -72,3 +83,16 @@ def _get_1d_funcs(self) -> Callable:
             raise ValueError(f"kernel embedding unknown.")
 
         return mean_func_1d
+    
+    def _get_1d_var_funcs(self) -> Callable:
+
+        var_func_1d_dict = {
+            "expquad-lebesgue": expquad_lebesgue_var_func_1d,
+        }
+
+        var_func_1d = var_func_1d_dict.get(self._kernel.name + "-" + self._measure.name, None)
+        if not var_func_1d:
+            pass
+            # raise ValueError(f"integrated kernel mean unknown.")
+
+        return var_func_1d
diff --git a/kernel_embedding_dictionary/embeddings/var_funcs_1d.py b/kernel_embedding_dictionary/embeddings/var_funcs_1d.py
new file mode 100644
index 0000000..9e1c7a2
--- /dev/null
+++ b/kernel_embedding_dictionary/embeddings/var_funcs_1d.py
@@ -0,0 +1,23 @@
+# Copyright 2025 The KED Authors. All Rights Reserved.
+# SPDX-License-Identifier: MIT
+
+
+import numpy as np
+from scipy.special import erf
+
+
+def expquad_lebesgue_var_func_1d(ell: float, lb: float, ub: float, density: float) -> np.ndarray:
+    """Compute the expected mean function for the exponential quadratic kernel with respect to the Lebesgue measure in 1D.
+    
+    Args:
+        ell: The length scale parameter.
+        lb: The lower bound of the integration interval.
+        ub: The upper bound of the integration interval.
+        density: The density of the Lebesgue measure.
+
+    Returns:
+        The expected mean function value.
+    """
+    exp_term = ell * np.sqrt(2/np.pi) * (np.exp(-0.5 * ((ub - lb) / ell) ** 2) - 1)
+    erf_term = (ub - lb) * (erf((ub - lb) / (ell * np.sqrt(2))))
+    return np.sqrt(2 * np.pi) * ell * density**2 * (exp_term + erf_term)
\ No newline at end of file
diff --git a/tests/kernel_embedding_dictionary/embeddings/test_expected_kmeans.py b/tests/kernel_embedding_dictionary/embeddings/test_expected_kmeans.py
new file mode 100644
index 0000000..f097a0a
--- /dev/null
+++ b/tests/kernel_embedding_dictionary/embeddings/test_expected_kmeans.py
@@ -0,0 +1,59 @@
+# Copyright 2025 The KED Authors. All Rights Reserved.
+# SPDX-License-Identifier: MIT
+
+import pytest
+import numpy as np
+
+from kernel_embedding_dictionary._get_embedding import get_embedding
+
+from scipy.integrate import quad
+
+
+def get_config_expquad_lebesgue_1d_standard():
+    ck = {"ndim": 1}
+    cm = {"ndim": 1}
+    return "expquad", "lebesgue", ck, cm
+
+
+def get_config_expquad_lebesgue_1d_values():
+    ck = {"ndim": 1, "lengthscales": [0.3]}
+    cm = {"ndim": 1, "bounds": [(-0.5, 2.5)], "normalize": True}  # test only works for normalized measures
+    return "expquad", "lebesgue", ck, cm
+
+
+@pytest.fixture()
+def config_expquad_lebesgue_1d_standard():
+    kn, mn, ck, cm = get_config_expquad_lebesgue_1d_standard()
+    ke = get_embedding(kernel_name=kn, measure_name=mn, kernel_config=ck, measure_config=cm)
+
+    def ke_mean_scalar(x):
+        return ke.mean(np.asarray(x).reshape(1, -1))
+
+    ekm, int_err = quad(ke_mean_scalar, 0, 1)
+    return kn, mn, ck, cm, ekm, int_err
+
+
+@pytest.fixture()
+def config_expquad_lebesgue_1d_values():
+    kn, mn, ck, cm = get_config_expquad_lebesgue_1d_values()
+    ke = get_embedding(kernel_name=kn, measure_name=mn, kernel_config=ck, measure_config=cm)
+
+    bounds = (cm["bounds"][0][0], cm["bounds"][0][1])
+
+    def ke_mean_scalar(x):
+        return ke.mean(np.asarray(x).reshape(1, -1)) / (bounds[1] - bounds[0])
+
+    ekm, int_err = quad(ke_mean_scalar, *bounds)
+    return kn, mn, ck, cm, ekm, int_err
+
+fixture_list = [
+    "config_expquad_lebesgue_1d_standard",
+    "config_expquad_lebesgue_1d_values",
+]
+
+@pytest.mark.parametrize("fixture_name", fixture_list)
+def test_expquad_lebesgue_mean_func_1d(fixture_name, request):
+    # Test cases for the expected mean function
+    kn, mn, ck, cm, num_ekm, int_err = request.getfixturevalue(fixture_name)
+    ke = get_embedding(kernel_name=kn, measure_name=mn, kernel_config=ck, measure_config=cm)
+    assert num_ekm == pytest.approx(ke.variance(), rel=int_err)
\ No newline at end of file