diff --git a/.test-conda-env-py3.yml b/.test-conda-env-py3.yml
index f9ede8fcf96f730ed9d14ce03d0504b81ab7ee13..eb7c3ddbc09c7ab0a82d08631f1562efadf53dc3 100644
--- a/.test-conda-env-py3.yml
+++ b/.test-conda-env-py3.yml
@@ -14,6 +14,7 @@ dependencies:
 - python=3
 - python-symengine=0.6.0
 - pyfmmlib
+- pyrsistent
 
 - pip
 - pip:
diff --git a/requirements.txt b/requirements.txt
index efe91f9487cb17f063b456973aabd909f153aae7..072cc9cef91b172b21b0a67484cab5274af04c3b 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,5 +1,6 @@
 numpy
 sympy==1.1.1
+pyrsistent
 git+https://github.com/inducer/pytools
 git+https://github.com/inducer/pymbolic
 git+https://github.com/inducer/islpy
diff --git a/sumpy/expansion/__init__.py b/sumpy/expansion/__init__.py
index fb19822b764fe21349f46aad3a40af6952599fd4..60647aac96ccebf6c896f0bfef8a33e5750e585b 100644
--- a/sumpy/expansion/__init__.py
+++ b/sumpy/expansion/__init__.py
@@ -354,7 +354,7 @@ class LinearPDEBasedExpansionTermsWrangler(ExpansionTermsWrangler):
         coeff_ident_enumerate_dict = {tuple(mi): i for
                                             (i, mi) in enumerate(mis)}
 
-        pde_dict = self.get_pde_as_diff_op().eq
+        pde_dict = self.get_pde_as_diff_op().mi_to_coeff
         for ident in pde_dict.keys():
             if ident not in coeff_ident_enumerate_dict:
                 # Order of the expansion is less than the order of the PDE.
diff --git a/sumpy/expansion/diff_op.py b/sumpy/expansion/diff_op.py
index 40dae7a434c0b14af8f1bc13a0908274f00d0b1b..73e85320e3835748112372ad259c0a7c8b50d3e3 100644
--- a/sumpy/expansion/diff_op.py
+++ b/sumpy/expansion/diff_op.py
@@ -22,6 +22,7 @@ OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
 THE SOFTWARE.
 """
 
+from pyrsistent import pmap
 from sumpy.tools import add_mi
 
 __doc__ = """
@@ -35,35 +36,35 @@ Differential operator interface
 class DifferentialOperator(object):
     r"""
     Represents a scalar, constant-coefficient DifferentialOperator of
-    dimension `dim`. It is represented by a dictionary. The dictionary
-    maps a multi-index given as a tuple to the coefficient.
+    dimension `dim`. It is represented by a frozen dictionary.
+    The dictionary maps a multi-index given as a tuple to the coefficient.
     This object is immutable.
     """
-    def __init__(self, dim, eq):
+    def __init__(self, dim, mi_to_coeff):
         """
         :arg dim: dimension of the DifferentialOperator
-        :arg eq: A dictionary mapping a multi-index to a coefficient
+        :arg mi_to_coeff: A dictionary mapping a multi-index to a coefficient
         """
         self.dim = dim
-        self.eq = eq
+        self.mi_to_coeff = mi_to_coeff
 
     def __mul__(self, param):
-        eq = {}
-        for k, v in self.eq.items():
-            eq[k] = v * param
-        return DifferentialOperator(self.dim, eq)
+        mi_to_coeff = {}
+        for k, v in self.mi_to_coeff.items():
+            mi_to_coeff[k] = v * param
+        return DifferentialOperator(self.dim, pmap(mi_to_coeff))
 
     __rmul__ = __mul__
 
     def __add__(self, other_pde):
         assert self.dim == other_pde.dim
-        res = self.eq.copy()
-        for k, v in other_pde.eq.items():
+        res = dict(self.mi_to_coeff)
+        for k, v in other_pde.mi_to_coeff.items():
             if k in res:
                 res[k] += v
             else:
                 res[k] = v
-        return DifferentialOperator(self.dim, res)
+        return DifferentialOperator(self.dim, pmap(res))
 
     __radd__ = __add__
 
@@ -71,12 +72,12 @@ class DifferentialOperator(object):
         return self + (-1)*other_pde
 
     def __repr__(self):
-        return f"DifferentialOperator({self.dim}, {repr(self.eq)})"
+        return f"DifferentialOperator({self.dim}, {repr(self.mi_to_coeff)})"
 
 
 def laplacian(pde):
     dim = pde.dim
-    res = DifferentialOperator(dim, {})
+    res = DifferentialOperator(dim, pmap())
     for j in range(dim):
         mi = [0]*dim
         mi[j] = 2
@@ -86,9 +87,9 @@ def laplacian(pde):
 
 def diff(pde, mi):
     res = {}
-    for eq_mi, v in pde.eq.items():
-        res[add_mi(eq_mi, mi)] = v
-    return DifferentialOperator(pde.dim, res)
+    for mi_to_coeff_mi, v in pde.mi_to_coeff.items():
+        res[add_mi(mi_to_coeff_mi, mi)] = v
+    return DifferentialOperator(pde.dim, pmap(res))
 
 
 def make_identity_diff_op(dim):
@@ -96,4 +97,4 @@ def make_identity_diff_op(dim):
     Returns the identity as a differential operator.
     """
     mi = tuple([0]*dim)
-    return DifferentialOperator(dim, {mi: 1})
+    return DifferentialOperator(dim, pmap({mi: 1}))