diff --git a/sumpy/expansion/local.py b/sumpy/expansion/local.py
index 483bc0eb16fcee7b43d54aedb4eaa8092f593620..809120776b79e518e4e8727dc9bcf94fa9082c55 100644
--- a/sumpy/expansion/local.py
+++ b/sumpy/expansion/local.py
@@ -133,7 +133,7 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase):
         bvec = bvec * rscale**-1
         result = sum(
                 coeff
-                * mi_power(bvec, mi)
+                * mi_power(bvec, mi, evaluate=False)
                 / mi_factorial(mi)
                 for coeff, mi in zip(
                         evaluated_coeffs, self.get_full_coefficient_identifiers()))
diff --git a/sumpy/symbolic.py b/sumpy/symbolic.py
index c7ed61e60c937cd18020b1ced0bd9629c5a0aa30..7a86958aebd7eea6c8054b11e61413314886ebe9 100644
--- a/sumpy/symbolic.py
+++ b/sumpy/symbolic.py
@@ -111,6 +111,14 @@ if not have_unevaluated_expr:
         return x
 
 
+if USE_SYMENGINE:
+    def unevaluated_pow(a, b):
+        return sym.Pow(a, b)
+else:
+    def unevaluated_pow(a, b):
+        return sym.Pow(a, b, evaluate=False)
+
+
 # {{{ debugging of sympy CSE via Maxima
 
 class _DerivativeKiller(IdentityMapperBase):
diff --git a/sumpy/tools.py b/sumpy/tools.py
index 469cb5e15238299092d55841f2bf43486bd11e0c..c826ef2fe8ea4d3f96e3b8b816dd3500f5efaae8 100644
--- a/sumpy/tools.py
+++ b/sumpy/tools.py
@@ -55,10 +55,15 @@ def mi_factorial(mi):
     return result
 
 
-def mi_power(vector, mi):
+def mi_power(vector, mi, evaluate=True):
     result = 1
     for mi_i, vec_i in zip(mi, vector):
-        result *= vec_i**mi_i
+        if mi_i == 1:
+            result *= vec_i
+        elif evaluate:
+            result *= vec_i**mi_i
+        else:
+            result *= sym.unevaluated_pow(vec_i, mi_i)
     return result