diff --git a/sumpy/fmm.py b/sumpy/fmm.py
index 519b12006f454f783056a4cb79b55d694eeb29d1..42564444afc604533f35dc268c8a2026c7fbc8db 100644
--- a/sumpy/fmm.py
+++ b/sumpy/fmm.py
@@ -323,8 +323,14 @@ class SumpyExpansionWrangler(ExpansionWranglerInterface):
         # YALE UNIV NEW HAVEN CT DEPT OF COMPUTER SCIENCE, 1988.
         # rscale that we use in sumpy is the inverse of the scaling used in the
         # paper and therefore we should use r / order. However empirically
-        # we have observed that 2r / order is better for numerical stability.
-        return r * 2 / order
+        # we have observed that 2r / order is better for numerical stability
+        # for Laplace and 4r / order for biharmonic kernel.
+        knl = self.tree_indep.get_base_kernel()
+        from sumpy.kernel import BiharmonicKernel
+        if isinstance(knl, BiharmonicKernel):
+            return r * 4 / order
+        else:
+            return r * 2 / order
 
     # {{{ data vector utilities
 
diff --git a/test/test_fmm.py b/test/test_fmm.py
index cbf66e6d125aa4f017551e931ad083cac7775276..4d6d0d480dbe1cc2fea190523bd291e421528b37 100644
--- a/test/test_fmm.py
+++ b/test/test_fmm.py
@@ -27,7 +27,8 @@ import numpy.linalg as la
 import pyopencl as cl
 from pyopencl.tools import (  # noqa
         pytest_generate_tests_for_pyopencl as pytest_generate_tests)
-from sumpy.kernel import LaplaceKernel, HelmholtzKernel, YukawaKernel
+from sumpy.kernel import (LaplaceKernel, HelmholtzKernel, YukawaKernel,
+    BiharmonicKernel)
 from sumpy.expansion.multipole import (
     VolumeTaylorMultipoleExpansion,
     H2DMultipoleExpansion, Y2DMultipoleExpansion,
@@ -234,6 +235,87 @@ def test_sumpy_fmm(ctx_factory, knl, local_expn_class, mpole_expn_class,
     pconv_verifier()
 
 
+@pytest.mark.parametrize("knl", [LaplaceKernel(2), BiharmonicKernel(2)])
+def test_coeff_magnitude_rscale(ctx_factory, knl):
+    """Checks that the rscale used keeps the coefficient magnitude
+    difference small
+    """
+    local_expn_class = LinearPDEConformingVolumeTaylorLocalExpansion
+    mpole_expn_class = LinearPDEConformingVolumeTaylorMultipoleExpansion
+
+    ctx = ctx_factory()
+    queue = cl.CommandQueue(ctx)
+
+    nsources = 1000
+    ntargets = 300
+    dtype = np.float64
+
+    from boxtree.tools import (
+            make_normal_particle_array as p_normal)
+
+    sources = p_normal(queue, nsources, knl.dim, dtype, seed=15)
+    offset = np.zeros(knl.dim)
+    offset[0] = 0.1
+
+    targets = (
+        p_normal(queue, ntargets, knl.dim, dtype, seed=18) + offset)
+
+    from boxtree import TreeBuilder
+    tb = TreeBuilder(ctx)
+
+    tree, _ = tb(queue, sources, targets=targets,
+            max_particles_in_box=30, debug=True)
+
+    from boxtree.traversal import FMMTraversalBuilder
+    tbuild = FMMTraversalBuilder(ctx)
+    trav, _ = tbuild(queue, tree, debug=True)
+
+    from pyopencl.clrandom import PhiloxGenerator
+    rng = PhiloxGenerator(ctx, seed=44)
+    weights = rng.uniform(queue, nsources, dtype=np.float64)
+
+    extra_kwargs = {}
+    dtype = np.float64
+    order = 10
+    if isinstance(knl, HelmholtzKernel):
+        extra_kwargs["k"] = 0.05
+        dtype = np.complex128
+
+    elif isinstance(knl, YukawaKernel):
+        extra_kwargs["lam"] = 2
+        dtype = np.complex128
+
+    from functools import partial
+    target_kernels = [knl]
+
+    tree_indep = SumpyTreeIndependentDataForWrangler(
+        ctx,
+        partial(mpole_expn_class, knl),
+        partial(local_expn_class, knl),
+        target_kernels)
+
+    def fmm_level_to_order(kernel, kernel_args, tree, lev):
+        return order
+
+    wrangler = SumpyExpansionWrangler(tree_indep, trav, dtype,
+        fmm_level_to_order=fmm_level_to_order,
+        kernel_extra_kwargs=extra_kwargs)
+
+    weights = wrangler.reorder_sources(weights)
+    (weights,) = wrangler.distribute_source_weights((weights,), None)
+
+    local_result, _ = wrangler.form_locals(
+        trav.level_start_target_or_target_parent_box_nrs,
+        trav.target_or_target_parent_boxes,
+        trav.from_sep_bigger_starts,
+        trav.from_sep_bigger_lists,
+        (weights,))
+
+    result = np.abs(wrangler.local_expansions_view(local_result, 5)[1][0])
+
+    assert np.max(result) / np.min(result) < 10**6
+
+
 def test_unified_single_and_double(ctx_factory):
     """
     Test that running one FMM for single layer + double layer gives the