diff --git a/examples/coeff-scaling-experiment.py b/examples/coeff-scaling-experiment.py
new file mode 100644
index 0000000000000000000000000000000000000000..3c14dd45be05cc3b3dfcbc7425f41ee4b7c0fc14
--- /dev/null
+++ b/examples/coeff-scaling-experiment.py
@@ -0,0 +1,55 @@
+import pyopencl as cl
+import sumpy.toys as t
+import numpy as np
+from sumpy.visualization import FieldPlotter
+import matplotlib.pyplot as plt
+
+
+# FIXME: Get rid of this once everything is working
+
+def main():
+    dim = 2
+
+    from sumpy.kernel import (  # noqa: F401
+            YukawaKernel, HelmholtzKernel, LaplaceKernel,
+            AxisTargetDerivative)
+    tctx = t.ToyContext(
+            cl.create_some_context(),
+            #LaplaceKernel(dim),
+            #AxisTargetDerivative(0, LaplaceKernel(dim)),
+            #YukawaKernel(dim), extra_source_kwargs={"lam": 5},
+            HelmholtzKernel(dim), extra_source_kwargs={"k": 0.3},
+            )
+
+    np.random.seed(12)
+    scale = 2**(-14)
+    #scale = 1
+    pts = np.random.rand(dim, 50) - 0.5
+    pt_src = t.PointSources(
+            tctx,
+            scale * pts,
+            np.ones(50))
+
+    mctr = scale*np.array([0., 0, 0])[:dim]
+    mexp = t.multipole_expand(pt_src, mctr, order=4, rscale=scale)
+
+    lctr = scale*np.array([2.5, 0, 0])[:dim]
+    #lexp = t.local_expand(mexp, lctr)
+    print(mexp.coeffs)
+    #print(lexp.coeffs)
+
+    diff = mexp - pt_src
+
+    diag = np.sqrt(dim)
+    print(t.l_inf(diff, scale*0.5*diag, center=lctr)
+            / t.l_inf(pt_src, scale*0.5*diag, center=lctr))
+    if 0:
+        fp = FieldPlotter(lexp.center, extent=scale*0.5)
+
+        t.logplot(fp, diff, cmap="jet", vmin=-16, vmax=-7)
+        plt.colorbar()
+        plt.show()
+
+
+if __name__ == "__main__":
+    main()
diff --git a/sumpy/e2p.py b/sumpy/e2p.py
index 6a77a5cc416af3c01df9786bac5fd236927e787c..ea5c382f66b9af73c51eaf995a3bc87e6d5fbeab 100644
--- a/sumpy/e2p.py
+++ b/sumpy/e2p.py
@@ -77,12 +77,15 @@ class E2PBase(KernelCacheWrapper):
         from sumpy.symbolic import make_sym_vector
         bvec = make_sym_vector("b", self.dim)
 
+        import sympy as sp
+        rscale = sp.Symbol("rscale")
+
         from sumpy.assignment_collection import SymbolicAssignmentCollection
         sac = SymbolicAssignmentCollection()
 
         coeff_exprs = [sym.Symbol("coeff%d" % i)
                 for i in range(len(self.expansion.get_coefficient_identifiers()))]
-        value = self.expansion.evaluate(coeff_exprs, bvec)
+        value = self.expansion.evaluate(coeff_exprs, bvec, rscale)
 
         result_names = [
             sac.assign_unique("result_%d_p" % i,
@@ -108,7 +111,8 @@ class E2PBase(KernelCacheWrapper):
         sympy_conv = SympyToPymbolicMapper()
         return [lp.Assignment(id=None,
                     assignee="kernel_scaling",
-                    expression=sympy_conv(self.expansion.kernel.get_scaling()),
+                    expression=sympy_conv(
+                        self.expansion.kernel.get_global_scaling_const()),
                     temp_var_type=lp.auto)]
 
     def get_cache_key(self):
@@ -165,6 +169,7 @@ class E2PFromSingleBox(E2PBase):
                     lp.GlobalArg("box_target_starts,box_target_counts_nonchild",
                         None, shape=None),
                     lp.GlobalArg("centers", None, shape="dim, naligned_boxes"),
+                    lp.ValueArg("rscale", None),
                     lp.GlobalArg("result", None, shape="nresults, ntargets",
                         dim_tags="sep,C"),
                     lp.GlobalArg("src_expansions", None,
@@ -205,7 +210,12 @@ class E2PFromSingleBox(E2PBase):
         """
         knl = self.get_cached_optimized_kernel()
 
-        return knl(queue, **kwargs)
+        centers = kwargs.pop("centers")
+        # "1" may be passed for rscale, which won't have its type
+        # meaningfully inferred. Make the type of rscale explicit.
+        rscale = centers.dtype.type(kwargs.pop("rscale"))
+
+        return knl(queue, centers=centers, rscale=rscale, **kwargs)
 
 # }}}
 
diff --git a/sumpy/expansion/__init__.py b/sumpy/expansion/__init__.py
index 4c76ea45d035e74df87b959606be4897d2deaa34..3810748220bd0b45d51f1e271a2bbe600dd5a373 100644
--- a/sumpy/expansion/__init__.py
+++ b/sumpy/expansion/__init__.py
@@ -494,8 +494,7 @@ class DefaultExpansionFactory(ExpansionFactoryBase):
                 and base_kernel.dim == 2):
             from sumpy.expansion.local import Y2DLocalExpansion
             return Y2DLocalExpansion
-        elif (isinstance(base_kernel.get_base_kernel(), HelmholtzKernel)
-                and base_kernel.dim == 3):
+        elif isinstance(base_kernel.get_base_kernel(), HelmholtzKernel):
             from sumpy.expansion.local import \
                     HelmholtzConformingVolumeTaylorLocalExpansion
             return HelmholtzConformingVolumeTaylorLocalExpansion
@@ -523,8 +522,7 @@ class DefaultExpansionFactory(ExpansionFactoryBase):
             from sumpy.expansion.multipole import (
                     LaplaceConformingVolumeTaylorMultipoleExpansion)
             return LaplaceConformingVolumeTaylorMultipoleExpansion
-        elif (isinstance(base_kernel.get_base_kernel(), HelmholtzKernel)
-                and base_kernel.dim == 3):
+        elif isinstance(base_kernel.get_base_kernel(), HelmholtzKernel):
             from sumpy.expansion.multipole import (
                     HelmholtzConformingVolumeTaylorMultipoleExpansion)
             return HelmholtzConformingVolumeTaylorMultipoleExpansion
diff --git a/sumpy/expansion/multipole.py b/sumpy/expansion/multipole.py
index 958bca6dbec6b8fdcbbb77545f039ae92a65d3d1..50d18142e27715837188d16d0a503a869f3440c6 100644
--- a/sumpy/expansion/multipole.py
+++ b/sumpy/expansion/multipole.py
@@ -25,6 +25,7 @@ THE SOFTWARE.
 from six.moves import range, zip
 import sumpy.symbolic as sym  # noqa
 
+from sumpy.tools import sympy_vec_subs
 from sumpy.expansion import (
     ExpansionBase, VolumeTaylorExpansion, LaplaceConformingVolumeTaylorExpansion,
     HelmholtzConformingVolumeTaylorExpansion)
@@ -53,12 +54,14 @@ class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase):
     Coefficients represent the terms in front of the kernel derivatives.
     """
 
-    def coefficients_from_source(self, avec, bvec):
+    def coefficients_from_source(self, avec, bvec, rscale):
         from sumpy.kernel import DirectionalSourceDerivative
         kernel = self.kernel
 
         from sumpy.tools import mi_power, mi_factorial
 
+        avec = avec/rscale
+
         if isinstance(kernel, DirectionalSourceDerivative):
             if kernel.get_base_kernel() is not kernel.inner_kernel:
                 raise NotImplementedError("more than one source derivative "
@@ -90,18 +93,24 @@ class VolumeTaylorMultipoleExpansionBase(MultipoleExpansionBase):
         return (
             self.derivative_wrangler.get_stored_mpole_coefficients_from_full(result))
 
-    def evaluate(self, coeffs, bvec):
-        taker = self.get_kernel_derivative_taker(bvec)
+    def evaluate(self, coeffs, bvec, rscale):
+        taker = self.get_kernel_derivative_taker(bvec, rscale)
         result = sym.Add(*tuple(
-                coeff * taker.diff(mi)
+                coeff
+                * self.kernel.adjust_proxy_expression(
+                    sympy_vec_subs(
+                        bvec, bvec/rscale,
+                        taker.diff(mi)),
+                    rscale,
+                    sum(mi))
+                * rscale**sum(mi)
                 for coeff, mi in zip(coeffs, self.get_coefficient_identifiers())))
-        return result
 
-    def get_kernel_derivative_taker(self, bvec):
-        ppkernel = self.kernel.postprocess_at_target(
-                self.kernel.get_expression(bvec), bvec)
+        return result
 
-        return self.derivative_wrangler.get_derivative_taker(ppkernel, bvec)
+    def get_kernel_derivative_taker(self, bvec, rscale):
+        return (self.derivative_wrangler.get_derivative_taker(
+            self.kernel.get_proxy_expression(bvec), bvec))
 
     def translate_from(self, src_expansion, src_coeff_exprs, dvec):
         if not isinstance(src_expansion, type(self)):
diff --git a/sumpy/kernel.py b/sumpy/kernel.py
index aed3d931d1fa18cb78825de0ad9ab5c9760c9ea7..f93a7575416daed6d3b48fd77e1fefb26a741945 100644
--- a/sumpy/kernel.py
+++ b/sumpy/kernel.py
@@ -31,8 +31,10 @@ from sumpy.symbolic import pymbolic_real_norm_2
 from pymbolic.primitives import make_sym_vector
 from pymbolic import var
 
-__doc__ = """
 
+# FIXME: Rename ProxyExpressionKernel back to ExpressionKernel once we're done.
+
+__doc__ = """
 Kernel interface
 ----------------
 
@@ -48,6 +50,7 @@ PDE kernels
 -----------
 
 .. autoclass:: LaplaceKernel
+.. autoclass:: BiharmonicKernel
 .. autoclass:: HelmholtzKernel
 .. autoclass:: YukawaKernel
 .. autoclass:: StokesletKernel
@@ -101,12 +104,20 @@ class Kernel(object):
     .. attribute:: is_complex_valued
     .. attribute:: dim
 
-        *dim* is allowed to be *None* if the dimensionality is not yet
-        known.
+    .. automethod:: get_base_kernel
+    .. automethod:: prepare_loopy_kernel
+    .. automethod:: get_code_transformer
+    .. automethod:: get_proxy_expression
+    .. automethod:: adjust_proxy_expression
+    .. automethod:: postprocess_at_source
+    .. automethod:: postprocess_at_target
+    .. automethod:: get_global_scaling_const
+    .. automethod:: get_args
+    .. automethod:: get_source_args
     """
 
     def __init__(self, dim=None):
-        self._dim = dim
+        self.dim = dim
 
     # {{{ hashing/pickling/equality
 
@@ -140,21 +151,10 @@ class Kernel(object):
         # Can't use trivial pickling: hash_value cache must stay unset
         assert len(self.init_arg_names) == len(state)
         for name, value in zip(self.init_arg_names, state):
-            if name == "dim":
-                name = "_dim"
-
             setattr(self, name, value)
 
     # }}}
 
-    @property
-    def dim(self):
-        if self._dim is None:
-            raise RuntimeError("the number of dimensions for this kernel "
-                    "has not yet been set")
-
-        return self._dim
-
     def get_base_kernel(self):
         """Return the kernel being wrapped by this one, or else
         *self*.
@@ -170,16 +170,62 @@ class Kernel(object):
     def get_code_transformer(self):
         """Return a function to postprocess the :mod:`pymbolic`
         expression generated from the result of
-        :meth:`get_expression` on the way to code generation.
+        :meth:`get_proxy_expression` on the way to code generation.
         """
         return lambda expr: expr
 
-    def get_expression(self, dist_vec):
-        """Return a :mod:`sympy` expression for the kernel.
-
-        :arg dist_vec: target - source
+    def get_proxy_expression(self, scaled_dist_vec):
+        r"""Return a proxy expression :math:`\tilde G(\vec \rho)` where
+        :math:`\vec \rho` is *scaled_dist_vec*. Here,
+        :math:`\vec \rho =\vec x/\alpha`,
+        where :math:`\alpha` is a scaling factor (often written as *rscale*)
+        used to ensure that :math:`\rho` has roughly unit magnitude and
+        :math:`\vec x=\text{target} - \text{source}`.  (In
+        an FMM, :math:`\alpha` would be chosen proportional to the
+        box size.) :math:`\tilde G` should be chosen so that (a)
+        no powers of :math:`\alpha` appear when symbolic derivatives
+        of it are taken with respect to :math:`\vec \rho` and
+        (b) the correct partial derivatives (with input scaling by :math:`\alpha`
+        can nonetheless be recovered) by :meth:`adjust_proxy_expression`.
+
+        Let :math:`\hat G` be the result of calling
+        :meth:`adjust_proxy_expression` with
+        :math:`\partial_{\vec \rho}^\mu \tilde G` and with
+        *nderivatives* set equal to :math:`|\mu|=\sum_i \mu_i`,
+        where :math:`\mu` is a multi-index indicating
+        how many derivatives are being taken along each input axis.
+
+        Then the following equality is required to hold:
+
+        .. math::
+
+            \hat G(\vec \rho) \alpha^{|\mu|}
+            = \partial_{\vec \rho}^\mu \tilde G(\vec \rho) \alpha^{|\mu|}
+            = \partial_{\vec x}^\mu G(\vec x).
+
+        In other words, :meth:`adjust_proxy_expression` is *not*
+        responsible for applying the scaling :math:`\alpha^{|\mu|}`
+        resulting from the chain rule (this is managed separately within
+        :mod:`sumpy` in an effort to control expansion coefficient scaling),
+        but it should compensate for any other effects of the
+        scaling of the input variable.
+
+        As a worked example, consider a kernel of :math:`\log(x)`. The
+        proxy expression would be :math:`\log(\rho)`.
+
+        For ``nderivatives == 0``, :meth:`adjust_proxy_expression`
+        would return :math:`\tilde G + \log(\alpha)`.
+
+        For ``nderivatives > 0``, :meth:`adjust_proxy_expression`
+        would return :math:`\tilde G/\alpha^{\text{nderivatives}}`.
+
+        In each case, the added terms merely compensate for the missing
+        input scaling.
+        """
+        raise NotImplementedError
 
-        (Assumes translation invariance of the kernel.)
+    def adjust_proxy_expression(self, expr, rscale, nderivatives):
+        r"""See :meth:`get_proxy_expression`.
         """
         raise NotImplementedError
 
@@ -187,6 +233,9 @@ class Kernel(object):
         """Transform a kernel evaluation or expansion expression in a place
         where the vector a (something - source) is known. ("something" may be
         an expansion center or a target.)
+
+        The typical use of this function is to apply source-variable
+        derivatives to the kernel.
         """
         return expr
 
@@ -194,11 +243,20 @@ class Kernel(object):
         """Transform a kernel evaluation or expansion expression in a place
         where the vector b (target - something) is known. ("something" may be
         an expansion center or a target.)
+
+        The typical use of this function is to apply target-variable
+        derivatives to the kernel.
         """
         return expr
 
-    def get_scaling(self):
-        """Return a global scaling of the kernel."""
+    def get_global_scaling_const(self):
+        """Return a global scaling constant of the kernel.
+         Typically, this ensures that the kernel is scaled so that
+         :math:`\mathcal L(G)(x)=C\delta(x)` with a constant of 1, where
+         :math:`\mathcal L` is the PDE operator associated with the kernel.
+         Not to be confused with *rscale*, which keeps expansion
+         coefficients benignly scaled.
+        """
         raise NotImplementedError
 
     def get_args(self):
@@ -218,32 +276,42 @@ class Kernel(object):
 
 
 class ExpressionKernel(Kernel):
+    def get_expression(self, dist_vec):
+        raise NotImplementedError
+
+
+class ProxyExpressionKernel(Kernel):
     is_complex_valued = False
 
-    init_arg_names = ("dim", "expression", "scaling", "is_complex_valued")
+    init_arg_names = ("dim", "expression", "kernel_scaling", "is_complex_valued")
 
-    def __init__(self, dim, expression, scaling, is_complex_valued):
+    def __init__(self, dim, proxy_expression, global_scaling_const,
+            is_complex_valued):
         """
         :arg expression: A :mod:`pymbolic` expression depending on
             variables *d_1* through *d_N* where *N* equals *dim*.
             (These variables match what is returned from
             :func:`pymbolic.primitives.make_sym_vector` with
             argument `"d"`.)
-        :arg scaling: A :mod:`pymbolic` expression for the scaling
-            of the kernel.
+        :arg global_scaling_const: A constant :mod:`pymbolic` expression for the
+            global scaling of the kernel. Typically, this ensures that
+            the kernel is scaled so that :math:`\mathcal L(G)(x)=C\delta(x)`
+            with a constant of 1, where :math:`\mathcal L` is the PDE
+            operator associated with the kernel. Not to be confused with
+            *rscale*, which keeps expansion coefficients benignly scaled.
         """
 
-        # expression and scaling are pymbolic objects because those pickle
-        # cleanly. D'oh, sympy!
+        # expression and global_scaling_const are pymbolic objects because
+        # those pickle cleanly. D'oh, sympy!
 
         Kernel.__init__(self, dim)
 
-        self.expression = expression
-        self.scaling = scaling
+        self.proxy_expression = proxy_expression
+        self.global_scaling_const = global_scaling_const
         self.is_complex_valued = is_complex_valued
 
     def __getinitargs__(self):
-        return (self._dim, self.expression, self.scaling,
+        return (self._dim, self.expression, self.global_scaling_const,
                 self.is_complex_valued)
 
     def __repr__(self):
@@ -252,36 +320,27 @@ class ExpressionKernel(Kernel):
         else:
             return "ExprKnl"
 
-    def get_expression(self, dist_vec):
-        if self.expression is None:
-            raise RuntimeError("expression in ExpressionKernel has not "
-                    "been determined yet (this could be due to a PDE kernel "
-                    "not having learned its dimensionality yet)")
-
+    def get_proxy_expression(self, scaled_dist_vec):
         from sumpy.symbolic import PymbolicToSympyMapperWithSymbols
-        expr = PymbolicToSympyMapperWithSymbols()(self.expression)
+        expr = PymbolicToSympyMapperWithSymbols()(self.proxy_expression)
 
-        if self.dim != len(dist_vec):
+        if self.dim != len(scaled_dist_vec):
             raise ValueError("dist_vec length does not match expected dimension")
 
         from sumpy.symbolic import Symbol
         expr = expr.subs(dict(
             (Symbol("d%d" % i), dist_vec_i)
-            for i, dist_vec_i in enumerate(dist_vec)
+            for i, dist_vec_i in enumerate(scaled_dist_vec)
             ))
 
         return expr
 
-    def get_scaling(self):
+    def get_global_scaling_const(self):
         """Return a global scaling of the kernel."""
 
-        if self.scaling is None:
-            raise RuntimeError("scaling in ExpressionKernel has not "
-                    "been determined yet (this could be due to a PDE kernel "
-                    "not having learned its dimensionality yet)")
-
         from sumpy.symbolic import PymbolicToSympyMapperWithSymbols
-        return PymbolicToSympyMapperWithSymbols()(self.scaling)
+        return PymbolicToSympyMapperWithSymbols()(
+                self.global_scaling_const)
 
     def update_persistent_hash(self, key_hash, key_builder):
         key_hash.update(type(self).__name__.encode("utf8"))
@@ -296,21 +355,21 @@ class ExpressionKernel(Kernel):
     mapper_method = "map_expression_kernel"
 
 
-one_kernel_2d = ExpressionKernel(
+one_kernel_2d = ProxyExpressionKernel(
         dim=2,
-        expression=1,
-        scaling=1,
+        proxy_expression=1,
+        global_scaling_const=1,
         is_complex_valued=False)
-one_kernel_3d = ExpressionKernel(
+one_kernel_3d = ProxyExpressionKernel(
         dim=3,
-        expression=1,
-        scaling=1,
+        proxy_expression=1,
+        global_scaling_const=1,
         is_complex_valued=False)
 
 
 # {{{ PDE kernels
 
-class LaplaceKernel(ExpressionKernel):
+class LaplaceKernel(ProxyExpressionKernel):
     init_arg_names = ("dim",)
 
     def __init__(self, dim=None):
@@ -324,15 +383,28 @@ class LaplaceKernel(ExpressionKernel):
             expr = 1/r
             scaling = 1/(4*var("pi"))
         else:
-            raise RuntimeError("unsupported dimensionality")
+            raise NotImplementedError("unsupported dimensionality")
 
-        ExpressionKernel.__init__(
-                self,
+        super(LaplaceKernel, self).__init__(
                 dim,
-                expression=expr,
-                scaling=scaling,
+                proxy_expression=expr,
+                global_scaling_const=scaling,
                 is_complex_valued=False)
 
+    def adjust_proxy_expression(self, expr, rscale, nderivatives):
+        if self.dim == 2:
+            if nderivatives == 0:
+                import sympy as sp
+                return expr + sp.log(rscale)
+            else:
+                return expr / rscale**nderivatives
+
+        elif self.dim == 3:
+            return expr / rscale**(nderivatives+1)
+
+        else:
+            raise NotImplementedError("unsupported dimensionality")
+
     def __getinitargs__(self):
         return (self._dim,)
 
@@ -354,6 +426,7 @@ class BiharmonicKernel(ExpressionKernel):
             expr = r**2 * var("log")(r)
             scaling = 1/(8*var("pi"))
         elif dim == 3:
+            raise RuntimeError("unsupported dimensionality")
             expr = r
             scaling = 1  # FIXME: Unknown
         else:
@@ -378,7 +451,7 @@ class BiharmonicKernel(ExpressionKernel):
     mapper_method = "map_biharmonic_kernel"
 
 
-class HelmholtzKernel(ExpressionKernel):
+class HelmholtzKernel(ProxyExpressionKernel):
     init_arg_names = ("dim", "helmholtz_k_name", "allow_evanescent")
 
     def __init__(self, dim=None, helmholtz_k_name="k",
@@ -403,16 +476,26 @@ class HelmholtzKernel(ExpressionKernel):
         else:
             raise RuntimeError("unsupported dimensionality")
 
-        ExpressionKernel.__init__(
-                self,
+        super(HelmholtzKernel, self).__init__(
                 dim,
-                expression=expr,
-                scaling=scaling,
+                proxy_expression=expr,
+                global_scaling_const=scaling,
                 is_complex_valued=True)
 
         self.helmholtz_k_name = helmholtz_k_name
         self.allow_evanescent = allow_evanescent
 
+    def adjust_proxy_expression(self, expr, rscale, nderivatives):
+        import sympy as sp
+        result = expr.subs(
+            sp.Symbol(self.helmholtz_k_name),
+            sp.Symbol(self.helmholtz_k_name) * rscale)
+
+        if self.dim == 2:
+            return result
+        else:
+            return result / rscale**(nderivatives+1)
+
     def __getinitargs__(self):
         return (self._dim, self.helmholtz_k_name,
                 self.allow_evanescent)
@@ -665,8 +748,11 @@ class KernelWrapper(Kernel):
     def is_complex_valued(self):
         return self.inner_kernel.is_complex_valued
 
-    def get_expression(self, dist_vec):
-        return self.inner_kernel.get_expression(dist_vec)
+    def get_proxy_expression(self, scaled_dist_vec):
+        return self.inner_kernel.get_proxy_expression(scaled_dist_vec)
+
+    def adjust_proxy_expression(self, expr, rscale, nderivatives):
+        return self.inner_kernel.adjust_proxy_expression(expr, rscale, nderivatives)
 
     def postprocess_at_source(self, expr, avec):
         return self.inner_kernel.postprocess_at_source(expr, avec)
@@ -674,8 +760,8 @@ class KernelWrapper(Kernel):
     def postprocess_at_target(self, expr, avec):
         return self.inner_kernel.postprocess_at_target(expr, avec)
 
-    def get_scaling(self):
-        return self.inner_kernel.get_scaling()
+    def get_global_scaling_const(self):
+        return self.inner_kernel.get_global_scaling_const()
 
     def get_code_transformer(self):
         return self.inner_kernel.get_code_transformer()
diff --git a/sumpy/p2e.py b/sumpy/p2e.py
index 6a8fcfa9d2ea8dd7a034d0adb1a7cc152c8db7bc..ce60eea961ae50be4dafe2f4c36cf678044ea7d9 100644
--- a/sumpy/p2e.py
+++ b/sumpy/p2e.py
@@ -73,13 +73,16 @@ class P2EBase(KernelCacheWrapper):
         from sumpy.symbolic import make_sym_vector
         avec = make_sym_vector("a", self.dim)
 
+        import sympy as sp
+        rscale = sp.Symbol("rscale")
+
         from sumpy.assignment_collection import SymbolicAssignmentCollection
         sac = SymbolicAssignmentCollection()
 
         coeff_names = [
                 sac.assign_unique("coeff%d" % i, coeff_i)
                 for i, coeff_i in enumerate(
-                    self.expansion.coefficients_from_source(avec, None))]
+                    self.expansion.coefficients_from_source(avec, None, rscale))]
 
         sac.run_global_cse()
 
@@ -93,7 +96,7 @@ class P2EBase(KernelCacheWrapper):
                 )
 
     def get_cache_key(self):
-        return (type(self).__name__, self.expansion)
+        return (type(self).__name__, self.name, self.expansion)
 
 # }}}
 
@@ -140,6 +143,7 @@ class P2EFromSingleBox(P2EBase):
                     lp.GlobalArg("box_source_starts,box_source_counts_nonchild",
                         None, shape=None),
                     lp.GlobalArg("centers", None, shape="dim, aligned_nboxes"),
+                    lp.ValueArg("rscale", None),
                     lp.GlobalArg("tgt_expansions", None,
                         shape=("nboxes", ncoeffs), offset=lp.auto),
                     lp.ValueArg("nboxes,aligned_nboxes,tgt_base_ibox", np.int32),
@@ -174,10 +178,16 @@ class P2EFromSingleBox(P2EBase):
         :arg centers:
         :arg sources:
         :arg strengths:
+        :arg rscale:
         """
+        centers = kwargs.pop("centers")
+        # "1" may be passed for rscale, which won't have its type
+        # meaningfully inferred. Make the type of rscale explicit.
+        rscale = centers.dtype.type(kwargs.pop("rscale"))
+
         knl = self.get_cached_optimized_kernel()
 
-        return knl(queue, **kwargs)
+        return knl(queue, centers=centers, rscale=rscale, **kwargs)
 
 # }}}
 
diff --git a/sumpy/p2p.py b/sumpy/p2p.py
index 9e07fbe01adffbe9a0732c46a80e87ee75b86e98..c7af82808ac4efe721ad4ddbf6e32662f13f538f 100644
--- a/sumpy/p2p.py
+++ b/sumpy/p2p.py
@@ -79,8 +79,12 @@ class P2PBase(KernelComputation, KernelCacheWrapper):
                 sac.assign_unique("knl%d" % i,
                     knl.postprocess_at_target(
                         knl.postprocess_at_source(
-                            knl.get_expression(dvec), dvec),
-                        dvec))
+                            knl.adjust_proxy_expression(
+                                knl.get_proxy_expression(dvec),
+                                rscale=1, nderivatives=0),
+                            dvec),
+                        dvec)
+                    )
                 for i, knl in enumerate(self.kernels)]
 
         sac.run_global_cse()
diff --git a/sumpy/tools.py b/sumpy/tools.py
index 2aa4732c02b6181d1324b9aa1dedee209cf949f9..6bda64f641f3b5ebd25683f38eb2df891b75cf04 100644
--- a/sumpy/tools.py
+++ b/sumpy/tools.py
@@ -271,7 +271,7 @@ class KernelComputation(object):
         return [
                 lp.Assignment(id=None,
                     assignee="knl_%d_scaling" % i,
-                    expression=sympy_conv(kernel.get_scaling()),
+                    expression=sympy_conv(kernel.get_global_scaling_const()),
                     temp_var_type=dtype)
                 for i, (kernel, dtype) in enumerate(
                     zip(self.kernels, self.value_dtypes))]
@@ -431,4 +431,14 @@ def my_syntactic_subs(expr, subst_dict):
         return expr
 
 
+def sympy_vec_subs(from_vec, to_vec, expr):
+    subs_pairs = []
+    assert (from_vec.rows, from_vec.cols) == (to_vec.rows, to_vec.cols)
+    for irow in range(from_vec.rows):
+        for icol in range(from_vec.cols):
+            subs_pairs.append((from_vec[irow, icol], to_vec[irow, icol]))
+
+    return expr.subs(subs_pairs)
+
+
 # vim: fdm=marker
diff --git a/sumpy/toys.py b/sumpy/toys.py
index e90bf5dce1b2f47cb1968f03bf11cccb33a2bf33..fae5eaf28f92b69eb25954e701c2495e21f7f4b5 100644
--- a/sumpy/toys.py
+++ b/sumpy/toys.py
@@ -127,7 +127,7 @@ class ToyContext(object):
 
 # {{{ helpers
 
-def _p2e(psource, center, order, p2e, expn_class, expn_kwargs):
+def _p2e(psource, center, rscale, order, p2e, expn_class, expn_kwargs):
     source_boxes = np.array([0], dtype=np.int32)
     box_source_starts = np.array([0], dtype=np.int32)
     box_source_counts_nonchild = np.array(
@@ -145,14 +145,15 @@ def _p2e(psource, center, order, p2e, expn_class, expn_kwargs):
             centers=centers,
             sources=psource.points,
             strengths=psource.weights,
+            rscale=rscale,
             nboxes=1,
             tgt_base_ibox=0,
 
             #flags="print_hl_cl",
             out_host=True, **toy_ctx.extra_source_kwargs)
 
-    return expn_class(toy_ctx, center, order, coeffs[0], derived_from=psource,
-            **expn_kwargs)
+    return expn_class(toy_ctx, center, rscale, order, coeffs[0],
+            derived_from=psource, **expn_kwargs)
 
 
 def _e2p(psource, targets, e2p):
@@ -176,6 +177,7 @@ def _e2p(psource, targets, e2p):
             box_target_starts=box_target_starts,
             box_target_counts_nonchild=box_target_counts_nonchild,
             centers=centers,
+            rscale=psource.rscale,
             targets=targets,
             #flags="print_hl_cl",
             out_host=True, **toy_ctx.extra_source_kwargs)
@@ -183,7 +185,7 @@ def _e2p(psource, targets, e2p):
     return pot
 
 
-def _e2e(psource, to_center, to_order, e2e, expn_class, expn_kwargs):
+def _e2e(psource, to_center, to_rscale, to_order, e2e, expn_class, expn_kwargs):
     toy_ctx = psource.toy_ctx
 
     target_boxes = np.array([1], dtype=np.int32)
@@ -217,7 +219,7 @@ def _e2e(psource, to_center, to_order, e2e, expn_class, expn_kwargs):
             #flags="print_hl_cl",
             out_host=True, **toy_ctx.extra_source_kwargs)
 
-    return expn_class(toy_ctx, to_center, to_order, to_coeffs[1],
+    return expn_class(toy_ctx, to_center, to_rscale, to_order, to_coeffs[1],
             derived_from=psource, **expn_kwargs)
 
 # }}}
@@ -319,10 +321,11 @@ class ExpansionPotentialSource(PotentialSource):
 
         Not used mathematically. Just for visualization, purely advisory.
     """
-    def __init__(self, toy_ctx, center, order, coeffs, derived_from,
+    def __init__(self, toy_ctx, center, rscale, order, coeffs, derived_from,
             radius=None, expn_style=None):
         super(ExpansionPotentialSource, self).__init__(toy_ctx)
         self.center = np.asarray(center)
+        self.rscale = rscale
         self.order = order
         self.coeffs = coeffs
 
@@ -380,12 +383,12 @@ class Product(PotentialExpressionNode):
 # }}}
 
 
-def multipole_expand(psource, center, order=None, **expn_kwargs):
+def multipole_expand(psource, center, order=None, rscale=1, **expn_kwargs):
     if isinstance(psource, PointSources):
         if order is None:
             raise ValueError("order may not be None")
 
-        return _p2e(psource, center, order, psource.toy_ctx.get_p2m(order),
+        return _p2e(psource, center, rscale, order, psource.toy_ctx.get_p2m(order),
                 MultipoleExpansion, expn_kwargs)
 
     elif isinstance(psource, MultipoleExpansion):
@@ -401,12 +404,12 @@ def multipole_expand(psource, center, order=None, **expn_kwargs):
                 % type(psource).__name__)
 
 
-def local_expand(psource, center, order=None, **expn_kwargs):
+def local_expand(psource, center, order=None, rscale=1, **expn_kwargs):
     if isinstance(psource, PointSources):
         if order is None:
             raise ValueError("order may not be None")
 
-        return _p2e(psource, center, order, psource.toy_ctx.get_p2l(order),
+        return _p2e(psource, center, rscale, order, psource.toy_ctx.get_p2l(order),
                 LocalExpansion, expn_kwargs)
 
     elif isinstance(psource, MultipoleExpansion):
diff --git a/sumpy/version.py b/sumpy/version.py
index 6333fc0cf50d9373f1beb568f773eba9c6422f25..72c4bd46b52b1ce24e238ca45c4241e3af8d7888 100644
--- a/sumpy/version.py
+++ b/sumpy/version.py
@@ -25,4 +25,4 @@ VERSION = (2016, 1)
 VERSION_STATUS = "beta1"
 VERSION_TEXT = ".".join(str(x) for x in VERSION) + VERSION_STATUS
 
-KERNEL_VERSION = 15
+KERNEL_VERSION = 16