diff --git a/pytential/symbolic/pde/scalar.py b/pytential/symbolic/pde/scalar.py
index 882b1b3c53cb53830bb94796b887a71494fdae32..0ddf9baff7620ba4fe62b3a58536257162cc5f84 100644
--- a/pytential/symbolic/pde/scalar.py
+++ b/pytential/symbolic/pde/scalar.py
@@ -26,6 +26,7 @@ __doc__ = """
 .. autoclass:: L2WeightedPDEOperator
 .. autoclass:: DirichletOperator
 .. autoclass:: NeumannOperator
+.. autoclass:: RobinOperator
 """
 
 
@@ -308,4 +309,117 @@ class NeumannOperator(L2WeightedPDEOperator):
 # }}}
 
 
+# {{{ robin
+
+class RobinOperator(L2WeightedPDEOperator):
+    """IE operator and field representation for solving Robin boundary
+    value problems with scalar kernels a * u + b * u_n = g (e.g.
+    :class:`sumpy.kernel.LaplaceKernel`)
+
+    .. automethod:: is_unique_only_up_to_constant
+    .. automethod:: representation
+    .. automethod:: operator
+    """
+
+    def __init__(self, kernel,
+            loc_sign, a, b,
+            use_l2_weighting=False, kernel_arguments=None):
+        """
+        :arg a: coefficient for the function value u.
+        :arg b: coefficient for the normal derivative u_n.
+        :arg loc_sign: +1 for exterior, -1 for interior
+        """
+
+        if not a**2 + b**2 > np.finfo(np.float64).eps:
+            raise ValueError("a and b cannot be both zero")
+
+        if a * b * loc_sign > np.finfo(np.float64).eps:
+            raise ValueError("illposed operator with a = %f, b = %f, loc_sign = %d"
+                    % (a, b, loc_sign))
+
+        assert loc_sign in [-1, 1]
+
+        from sumpy.kernel import LaplaceKernel
+        if not isinstance(kernel, LaplaceKernel):
+            raise NotImplementedError("Unsupported kernel %s" % str(kernel))
+
+        if kernel_arguments is None:
+            kernel_arguments = {}
+        self.kernel_arguments = kernel_arguments
+
+        self.a = a
+        self.b = b
+        self.alpha = b / (a**2 + b**2)
+        self.beta = -1 * a / (a**2 + b**2)
+        self.loc_sign = loc_sign
+
+        L2WeightedPDEOperator.__init__(self, kernel, use_l2_weighting)
+
+    def is_unique_only_up_to_constant(self):
+        from sumpy.kernel import LaplaceKernel
+        laplace_interior_neumann = (
+            isinstance(self.kernel, LaplaceKernel)
+            and self.loc_sign < 0
+            and self.a == 0)
+        laplace_exterior_dirichlet = (
+            isinstance(self.kernel, LaplaceKernel)
+            and self.loc_sign > 0
+            and self.b == 0)
+        return laplace_interior_neumann or laplace_exterior_dirichlet
+
+    def representation(self, u, map_potentials=None, qbx_forced_limit=None,
+            **kwargs):
+        sqrt_w = self.get_sqrt_weight()
+        inv_sqrt_w_u = cse(u/sqrt_w)
+
+        if map_potentials is None:
+            def map_potentials(x):  # pylint:disable=function-redefined
+                return x
+
+        kwargs["qbx_forced_limit"] = qbx_forced_limit
+        kwargs["kernel_arguments"] = self.kernel_arguments
+
+        return (
+                self.alpha * map_potentials(
+                    sym.S(self.kernel, inv_sqrt_w_u, **kwargs))
+                + self.beta * map_potentials(
+                    sym.D(self.kernel, inv_sqrt_w_u, **kwargs))
+                )
+
+    def operator(self, u):
+        sqrt_w = self.get_sqrt_weight()
+        inv_sqrt_w_u = cse(u/sqrt_w)
+
+        knl = self.kernel
+
+        knl_kwargs = {}
+        knl_kwargs["kernel_arguments"] = self.kernel_arguments
+
+        if self.is_unique_only_up_to_constant():
+            # The interior Neumann operator in this representation
+            # has a nullspace. The mean of the density must be matched
+            # to the desired solution separately. As is, this operator
+            # returns a mean that is not well-specified.
+
+            amb_dim = self.kernel.dim
+            ones_contribution = (
+                    sym.Ones() * sym.mean(amb_dim, amb_dim-1, inv_sqrt_w_u))
+        else:
+            ones_contribution = 0
+
+        return (-self.loc_sign * 0.5 * u
+                + sqrt_w * (
+                    self.a * self.alpha * sym.S(
+                        knl, inv_sqrt_w_u, qbx_forced_limit="avg", **knl_kwargs)
+                    + self.a * self.beta * sym.D(
+                        knl, inv_sqrt_w_u, qbx_forced_limit="avg", **knl_kwargs)
+                    + self.b * self.alpha * sym.Sp(
+                        knl, inv_sqrt_w_u, qbx_forced_limit="avg", **knl_kwargs)
+                    + self.b * self.beta * sym.Dp(
+                        knl, inv_sqrt_w_u, qbx_forced_limit="avg", **knl_kwargs)
+                    + ones_contribution
+                    ))
+
+# }}}
+
 # vim: foldmethod=marker