diff --git a/chapman_enskog/chapman_enskog.py b/chapman_enskog/chapman_enskog.py
index 33b82e1bd5c8ec17ccc1f64543224ae36437679b..b0ef2b38abf3bd1b74272c01f180289fe792f53c 100644
--- a/chapman_enskog/chapman_enskog.py
+++ b/chapman_enskog/chapman_enskog.py
@@ -367,11 +367,12 @@ def getTaylorExpandedLbEquation(pdfSymbolName="f", pdfsAfterCollisionOperator="\
 
     functions = [pdf, collidedPdf]
     eq_4_5 = taylorOperator - dt * collidedPdf
-    applied_eq_4_5 = expandUsingLinearity(DiffOperator.apply(eq_4_5, pdf), functions)
+    applied_eq_4_5 = expandUsingLinearity(DiffOperator.apply(eq_4_5, pdf, applyToConstants=False), functions)
 
     if shift:
         operator = ((dt / 2) * (Dt + c.dot(Dx))).expand()
-        opTimesEq_4_5 = expandUsingLinearity(DiffOperator.apply(operator, applied_eq_4_5), functions).expand()
+        opTimesEq_4_5 = expandUsingLinearity(DiffOperator.apply(operator, applied_eq_4_5, applyToConstants=False),
+                                             functions).expand()
         opTimesEq_4_5 = normalizeDiffOrder(opTimesEq_4_5, functions)
         eq_4_7 = (applied_eq_4_5 - opTimesEq_4_5).subs(dt ** (taylorOrder+1), 0)
     else:
@@ -671,7 +672,7 @@ class SteadyStateChapmanEnskogAnalysis(object):
         c = sp.Matrix([expandedSymbol("c", subscript=i) for i in range(dim)])
         Dx = sp.Matrix([DiffOperator(target=l) for l in range(dim)])
         differentialOperator = sum((dt * eps * c.dot(Dx)) ** n / sp.factorial(n) for n in range(1, order + 1))
-        taylorExpansion = DiffOperator.apply(differentialOperator.expand(), f)
+        taylorExpansion = DiffOperator.apply(differentialOperator.expand(), f, applyToConstants=False)
         epsDict = useChapmanEnskogAnsatz(taylorExpansion,
                                          spatialDerivativeOrders=None,  # do not expand the differential operator itself
                                          pdfs=(['f', 0, order + 1],))  # expand only the 'f' terms
diff --git a/phasefield/analytical.py b/phasefield/analytical.py
index 977096920677bb822022540d73e963c8ceac70d8..95c5bff24d6325e51ebe02514005d0636523679c 100644
--- a/phasefield/analytical.py
+++ b/phasefield/analytical.py
@@ -195,7 +195,8 @@ def chemicalPotentialsFromFreeEnergy(freeEnergy, orderParameters=None):
         orderParameters.sort(key=lambda e: e.name)
         orderParameters = orderParameters[:-1]
     constants = [s for s in syms if s not in orderParameters]
-    return sp.Matrix([functionalDerivative(freeEnergy, op, constants) for op in orderParameters])
+    return sp.Matrix([expandUsingLinearity(functionalDerivative(freeEnergy, op),constants=constants)
+                      for op in orderParameters])
 
 
 def forceFromPhiAndMu(orderParameters, dim, mu=None):