diff --git a/sumpy/expansion/__init__.py b/sumpy/expansion/__init__.py
index 7ab2125c7cf091f52a32fc4a168b98302d063710..22223b47fca42e15d3c3768216d8566d36bd9146 100644
--- a/sumpy/expansion/__init__.py
+++ b/sumpy/expansion/__init__.py
@@ -28,13 +28,12 @@ import numpy as np
 import logging
 from pytools import memoize_method
 import sumpy.symbolic as sym
-from sumpy.tools import MiDerivativeTaker
 from collections import defaultdict
 
 
 __doc__ = """
 .. autoclass:: ExpansionBase
-.. autoclass:: LinearRecurrenceBasedDerivativeWrangler
+.. autoclass:: LinearRecurrenceBasedExpansionTermsWrangler
 
 Expansion Factories
 ^^^^^^^^^^^^^^^^^^^
@@ -159,9 +158,9 @@ class ExpansionBase(object):
 # }}}
 
 
-# {{{ derivative wrangler
+# {{{ expansion terms wrangler
 
-class DerivativeWrangler(object):
+class ExpansionTermsWrangler(object):
 
     def __init__(self, order, dim):
         self.order = order
@@ -178,9 +177,6 @@ class DerivativeWrangler(object):
             rscale):
         raise NotImplementedError
 
-    def get_derivative_taker(self, expr, var_list):
-        raise NotImplementedError
-
     @memoize_method
     def get_full_coefficient_identifiers(self):
         """
@@ -193,18 +189,21 @@ class DerivativeWrangler(object):
         res = sorted(gnitstam(self.order, self.dim), key=sum)
         return res
 
-    @memoize_method
-    def copy(self, order):
-        raise NotImplementedError
+    def copy(self, **kwargs):
+        order = kwargs.pop('order', self.order)
+        dim = kwargs.pop('dim', self.dim)
+
+        if kwargs:
+            raise TypeError("unexpected keyword arguments '%s'"
+                % ", ".join(kwargs))
 
+        return type(self)(order, dim)
 
-class FullDerivativeWrangler(DerivativeWrangler):
 
-    def get_derivative_taker(self, expr, dvec):
-        return MiDerivativeTaker(expr, dvec)
+class FullExpansionTermsWrangler(ExpansionTermsWrangler):
 
     get_coefficient_identifiers = (
-            DerivativeWrangler.get_full_coefficient_identifiers)
+            ExpansionTermsWrangler.get_full_coefficient_identifiers)
 
     def get_full_kernel_derivatives_from_stored(self, stored_kernel_derivatives,
             rscale):
@@ -213,17 +212,6 @@ class FullDerivativeWrangler(DerivativeWrangler):
     get_stored_mpole_coefficients_from_full = (
             get_full_kernel_derivatives_from_stored)
 
-    @memoize_method
-    def copy(self, **kwargs):
-        order = kwargs.pop('order', self.order)
-        dim = kwargs.pop('dim', self.dim)
-
-        if kwargs:
-            raise TypeError("unexpected keyword arguments '%s'"
-                % ", ".join(kwargs))
-
-        return type(self)(order, dim)
-
 
 # {{{ sparse matrix-vector multiplication
 
@@ -260,7 +248,7 @@ def _spmv(spmat, x, sparse_vectors):
 # }}}
 
 
-class LinearRecurrenceBasedDerivativeWrangler(DerivativeWrangler):
+class LinearRecurrenceBasedExpansionTermsWrangler(ExpansionTermsWrangler):
     """
     .. automethod:: __init__
     .. automethod:: get_pde_dict
@@ -271,13 +259,14 @@ class LinearRecurrenceBasedDerivativeWrangler(DerivativeWrangler):
         r"""
         :param order: order of the expansion
         :param dim: number of dimensions
-        :param deriv_multiplier: a symbolic expression that is used to 'normalize out'
-            constant coefficients in the PDE in 
-            :func:`~LinearRecurrenceBasedDerivativeWrangler.get_pde_dict`, so that the
-            Taylor coefficient with multi-index :math:`\nu` as seen by that representation
-            of the PDE is :math:`\text{coeff} / {\text{deriv\_multiplier}^{|\nu|}}`.
+        :param deriv_multiplier: a symbolic expression that is used to
+            'normalize out' constant coefficients in the PDE in
+            :func:`~LinearRecurrenceBasedExpansionTermsWrangler.get_pde_dict`, so that
+            the Taylor coefficient with multi-index :math:`\nu` as seen by that
+            representation of the PDE is :math:`\text{coeff} /
+            {\text{deriv\_multiplier}^{|\nu|}}`.
         """
-        super(LinearRecurrenceBasedDerivativeWrangler, self).__init__(order, dim)
+        super(LinearRecurrenceBasedExpansionTermsWrangler, self).__init__(order, dim)
         self.deriv_multiplier = deriv_multiplier
 
     def get_coefficient_identifiers(self):
@@ -451,15 +440,11 @@ class LinearRecurrenceBasedDerivativeWrangler(DerivativeWrangler):
         """
         raise NotImplementedError
 
-    def get_derivative_taker(self, expr, var_list):
-        from sumpy.tools import MiDerivativeTaker
-        return MiDerivativeTaker(expr, var_list)
-
 
-class LaplaceDerivativeWrangler(LinearRecurrenceBasedDerivativeWrangler):
+class LaplaceExpansionTermsWrangler(LinearRecurrenceBasedExpansionTermsWrangler):
 
     def __init__(self, order, dim):
-        super(LaplaceDerivativeWrangler, self).__init__(order, dim, 1)
+        super(LaplaceExpansionTermsWrangler, self).__init__(order, dim, 1)
 
     def get_pde_dict(self):
         pde_dict = {}
@@ -479,11 +464,11 @@ class LaplaceDerivativeWrangler(LinearRecurrenceBasedDerivativeWrangler):
         return idx
 
 
-class HelmholtzDerivativeWrangler(LinearRecurrenceBasedDerivativeWrangler):
+class HelmholtzExpansionTermsWrangler(LinearRecurrenceBasedExpansionTermsWrangler):
 
     def __init__(self, order, dim, helmholtz_k_name):
         multiplier = sym.Symbol(helmholtz_k_name)
-        super(HelmholtzDerivativeWrangler, self).__init__(order, dim, multiplier)
+        super(HelmholtzExpansionTermsWrangler, self).__init__(order, dim, multiplier)
 
     def get_pde_dict(self, **kwargs):
         pde_dict = {}
@@ -511,31 +496,33 @@ class HelmholtzDerivativeWrangler(LinearRecurrenceBasedDerivativeWrangler):
 class VolumeTaylorExpansionBase(object):
 
     @classmethod
-    def get_or_make_derivative_wrangler(cls, *key):
+    def get_or_make_expansion_terms_wrangler(cls, *key):
         """
-        This stores the derivative wrangler at the class attribute level because
-        precomputing the derivative wrangler may become expensive.
+        This stores the expansion terms wrangler at the class attribute level
+        because recreating the expansion terms wrangler implicitly empties its
+        caches.
         """
         try:
-            wrangler = cls.derivative_wrangler_cache[key]
+            wrangler = cls.expansion_terms_wrangler_cache[key]
         except KeyError:
-            wrangler = cls.derivative_wrangler_class(*key)
-            cls.derivative_wrangler_cache[key] = wrangler
+            wrangler = cls.expansion_terms_wrangler_class(*key)
+            cls.expansion_terms_wrangler_cache[key] = wrangler
 
         return wrangler
 
     @property
-    def derivative_wrangler(self):
-        return self.get_or_make_derivative_wrangler(*self.derivative_wrangler_key)
+    def expansion_terms_wrangler(self):
+        return self.get_or_make_expansion_terms_wrangler(
+                *self.expansion_terms_wrangler_key)
 
     def get_coefficient_identifiers(self):
         """
         Returns the identifiers of the coefficients that actually get stored.
         """
-        return self.derivative_wrangler.get_coefficient_identifiers()
+        return self.expansion_terms_wrangler.get_coefficient_identifiers()
 
     def get_full_coefficient_identifiers(self):
-        return self.derivative_wrangler.get_full_coefficient_identifiers()
+        return self.expansion_terms_wrangler.get_full_coefficient_identifiers()
 
     @property
     @memoize_method
@@ -549,33 +536,33 @@ class VolumeTaylorExpansionBase(object):
 
 class VolumeTaylorExpansion(VolumeTaylorExpansionBase):
 
-    derivative_wrangler_class = FullDerivativeWrangler
-    derivative_wrangler_cache = {}
+    expansion_terms_wrangler_class = FullExpansionTermsWrangler
+    expansion_terms_wrangler_cache = {}
 
     # not user-facing, be strict about having to pass use_rscale
     def __init__(self, kernel, order, use_rscale):
-        self.derivative_wrangler_key = (order, kernel.dim)
+        self.expansion_terms_wrangler_key = (order, kernel.dim)
 
 
 class LaplaceConformingVolumeTaylorExpansion(VolumeTaylorExpansionBase):
 
-    derivative_wrangler_class = LaplaceDerivativeWrangler
-    derivative_wrangler_cache = {}
+    expansion_terms_wrangler_class = LaplaceExpansionTermsWrangler
+    expansion_terms_wrangler_cache = {}
 
     # not user-facing, be strict about having to pass use_rscale
     def __init__(self, kernel, order, use_rscale):
-        self.derivative_wrangler_key = (order, kernel.dim)
+        self.expansion_terms_wrangler_key = (order, kernel.dim)
 
 
 class HelmholtzConformingVolumeTaylorExpansion(VolumeTaylorExpansionBase):
 
-    derivative_wrangler_class = HelmholtzDerivativeWrangler
-    derivative_wrangler_cache = {}
+    expansion_terms_wrangler_class = HelmholtzExpansionTermsWrangler
+    expansion_terms_wrangler_cache = {}
 
     # not user-facing, be strict about having to pass use_rscale
     def __init__(self, kernel, order, use_rscale):
         helmholtz_k_name = kernel.get_base_kernel().helmholtz_k_name
-        self.derivative_wrangler_key = (order, kernel.dim, helmholtz_k_name)
+        self.expansion_terms_wrangler_key = (order, kernel.dim, helmholtz_k_name)
 
 # }}}
 
diff --git a/sumpy/expansion/local.py b/sumpy/expansion/local.py
index 0336ef43fc5f754782c632376605f7d9d54cc839..716fa93f4cedc71703f86e325421d69541e1df31 100644
--- a/sumpy/expansion/local.py
+++ b/sumpy/expansion/local.py
@@ -128,7 +128,7 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase):
     def evaluate(self, coeffs, bvec, rscale):
         from sumpy.tools import mi_power, mi_factorial
         evaluated_coeffs = (
-            self.derivative_wrangler.get_full_kernel_derivatives_from_stored(
+            self.expansion_terms_wrangler.get_full_kernel_derivatives_from_stored(
                 coeffs, rscale))
         bvec = bvec * rscale**-1
         result = sum(
@@ -167,50 +167,52 @@ class VolumeTaylorLocalExpansionBase(LocalExpansionBase):
 
             from sumpy.tools import add_mi
 
-            src_max_sum = max(sum(mi) for mi in
-                              src_expansion.get_coefficient_identifiers())
-            assert src_max_sum == src_expansion.order
-            tgt_max_sum = max(sum(mi) for mi in
-                              self.get_coefficient_identifiers())
-            assert tgt_max_sum == tgt_expansion.order            
+            # Create a expansion terms wrangler for derivatives up to order
+            # (tgt order)+(src order).
+            tgtplusderiv_exp_terms_wrangler = \
+                src_expansion.expansion_terms_wrangler.copy(
+                        order=self.order + src_expansion.order)
+            tgtplusderiv_coeff_ids = \
+                tgtplusderiv_exp_terms_wrangler.get_coefficient_identifiers()
+            tgtplusderiv_full_coeff_ids = \
+                tgtplusderiv_exp_terms_wrangler.get_full_coefficient_identifiers()
 
-            max_sum = src_max_sum + tgt_max_sum
-            new_deriv_wrangler = \
-                src_expansion.derivative_wrangler.copy(order=max_sum)
-            new_coeffs = new_deriv_wrangler.get_coefficient_identifiers()
-            new_full_coeffs = new_deriv_wrangler.get_full_coefficient_identifiers()
-
-            ident_to_index = dict((ident, i) for i, ident in
-                                enumerate(new_full_coeffs))
+            tgtplusderiv_ident_to_index = dict((ident, i) for i, ident in
+                                enumerate(tgtplusderiv_full_coeff_ids))
 
             result = []
-            for deriv in self.get_coefficient_identifiers():
-                local_result = []
+            for lexp_mi in self.get_coefficient_identifiers():
+                lexp_mi_terms = []
 
-                full_coeffs = [0] * len(new_full_coeffs)
+                # Embed the source coefficients in the coefficient array
+                # for the (tgt order)+(src order) wrangler, offset by lexp_mi.
+                embedded_coeffs = [0] * len(tgtplusderiv_full_coeff_ids)
                 for coeff, term in zip(
                         src_coeff_exprs,
                         src_expansion.get_coefficient_identifiers()):
-                    full_coeffs[ident_to_index[add_mi(deriv, term)]] = coeff
+                    embedded_coeffs[
+                            tgtplusderiv_ident_to_index[add_mi(lexp_mi, term)]] \
+                                    = coeff
 
-                stored_coeffs = \
-                    new_deriv_wrangler.get_stored_mpole_coefficients_from_full(
-                        full_coeffs, src_rscale)
+                stored_coeffs = tgtplusderiv_exp_terms_wrangler \
+                        .get_stored_mpole_coefficients_from_full(
+                                embedded_coeffs, src_rscale)
 
                 for i, coeff in enumerate(stored_coeffs):
                     if coeff == 0:
                         continue
-                    nderivatives_for_scaling = sum(new_coeffs[i])-sum(deriv)
+                    nderivatives_for_scaling = \
+                            sum(tgtplusderiv_coeff_ids[i])-sum(lexp_mi)
                     kernel_deriv = (
                             src_expansion.get_scaled_multipole(
-                                taker.diff(new_coeffs[i]),
+                                taker.diff(tgtplusderiv_coeff_ids[i]),
                                 dvec, src_rscale,
-                                nderivatives=sum(new_coeffs[i]),
+                                nderivatives=sum(tgtplusderiv_coeff_ids[i]),
                                 nderivatives_for_scaling=nderivatives_for_scaling))
 
-                    local_result.append(
-                            coeff * kernel_deriv * tgt_rscale**sum(deriv))
-                result.append(sym.Add(*local_result))
+                    lexp_mi_terms.append(
+                            coeff * kernel_deriv * tgt_rscale**sum(lexp_mi))
+                result.append(sym.Add(*lexp_mi_terms))
 
         else:
             from sumpy.tools import MiDerivativeTaker
diff --git a/sumpy/expansion/multipole.py b/sumpy/expansion/multipole.py
index 92ae80f3308b4bf8820fbc00da9229725058e3e1..25488952e4b7425371ecf23aefa261c0c70900f5 100644
--- a/sumpy/expansion/multipole.py
+++ b/sumpy/expansion/multipole.py
@@ -94,7 +94,7 @@ class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase):
                     mi_power(avec, mi) / mi_factorial(mi)
                     for mi in self.get_full_coefficient_identifiers()]
         return (
-            self.derivative_wrangler.get_stored_mpole_coefficients_from_full(
+            self.expansion_terms_wrangler.get_stored_mpole_coefficients_from_full(
                 result, rscale))
 
     def get_scaled_multipole(self, expr, bvec, rscale, nderivatives,
@@ -127,8 +127,8 @@ class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase):
         return result
 
     def get_kernel_derivative_taker(self, bvec):
-        return (self.derivative_wrangler.get_derivative_taker(
-            self.kernel.get_expression(bvec), bvec))
+        from sumpy.tools import MiDerivativeTaker
+        return MiDerivativeTaker(self.kernel.get_expression(bvec), bvec)
 
     def translate_from(self, src_expansion, src_coeff_exprs, src_rscale,
             dvec, tgt_rscale):
@@ -188,7 +188,7 @@ class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase):
 
         logger.info("building translation operator: done")
         return (
-            self.derivative_wrangler.get_stored_mpole_coefficients_from_full(
+            self.expansion_terms_wrangler.get_stored_mpole_coefficients_from_full(
                 result, tgt_rscale))