diff --git a/lumicks/pylake/fitting/parameters.py b/lumicks/pylake/fitting/parameters.py
index 18f841315..c9352d9fd 100644
--- a/lumicks/pylake/fitting/parameters.py
+++ b/lumicks/pylake/fitting/parameters.py
@@ -43,6 +43,7 @@ def __init__(
         fixed=False,
         shared=False,
         unit=None,
+        stderr=None,
     ):
         """Model parameter
 
@@ -90,7 +91,7 @@ def __init__(
         from the data. See also: :meth:`~lumicks.pylake.FdFit.profile_likelihood()`.
         """
 
-        self.stderr = None
+        self.stderr = stderr
         """Standard error of this parameter.
 
         Standard errors are calculated after fitting the model. These asymptotic errors are based
diff --git a/lumicks/pylake/fitting/profile_likelihood.py b/lumicks/pylake/fitting/profile_likelihood.py
index db973e1e3..715dfb534 100644
--- a/lumicks/pylake/fitting/profile_likelihood.py
+++ b/lumicks/pylake/fitting/profile_likelihood.py
@@ -624,7 +624,7 @@ def chi2(self):
     def p(self):
         return self.parameters[:, self.profile_info.profiled_parameter_index]
 
-    def plot(self, *, significance_level=None, **kwargs):
+    def plot(self, *, significance_level=None, std_err=None, **kwargs):
         """Plot profile likelihood
 
         Parameters
@@ -632,9 +632,19 @@ def plot(self, *, significance_level=None, **kwargs):
         significance_level : float, optional
             Desired significance level  (resulting in a 100 * (1 - alpha)% confidence interval) to
             plot. Default is the significance level specified when the profile was generated.
+        std_err : float | None
+            If provided, also make a quadratic plot based on a standard error.
         """
         import matplotlib.pyplot as plt
 
+        if std_err:
+            x = np.arange(-3 * std_err, 3 * std_err, 0.1 * std_err)
+            plt.plot(
+                self.p[np.argmin(self.chi2)] + x,
+                self.profile_info.minimum_chi2 + x**2 / (2 * std_err**2),
+                "k--",
+            )
+
         dash_length = 5
         plt.plot(self.p, self.chi2, **kwargs)
 
diff --git a/lumicks/pylake/population/dwelltime.py b/lumicks/pylake/population/dwelltime.py
index 35f16c59d..0644a10b8 100644
--- a/lumicks/pylake/population/dwelltime.py
+++ b/lumicks/pylake/population/dwelltime.py
@@ -99,8 +99,10 @@ def fit_func(params, lb, ub, fitted):
         )
         parameters = Params(
             **{
-                key: Parameter(param, lower_bound=lb, upper_bound=ub)
-                for key, param, (lb, ub) in zip(keys, dwelltime_model._parameters, bounds)
+                key: Parameter(param, lower_bound=lb, upper_bound=ub, stderr=std_err)
+                for key, param, (lb, ub), std_err in zip(
+                    keys, dwelltime_model._parameters, bounds, dwelltime_model._std_errs
+                )
             }
         )
 
@@ -195,7 +197,7 @@ def n_components(self):
         """Number of components in the model."""
         return self.model.n_components
 
-    def plot(self, alpha=None):
+    def plot(self, alpha=None, *, with_stderr=False, **kwargs):
         """Plot the profile likelihoods for the parameters of a model.
 
         Confidence interval is indicated by the region where the profile crosses the chi squared
@@ -207,16 +209,26 @@ def plot(self, alpha=None):
             Significance level. Confidence intervals are calculated as 100 * (1 - alpha)%. The
             default value of `None` results in using the significance level applied when
             profiling (default: 0.05).
+        with_stderr : bool
+            Also show bounds based on standard errors.
         """
         import matplotlib.pyplot as plt
 
+        std_errs = self.model._std_errs[~np.isnan(self.model._std_errs)]
         if self.n_components == 1:
-            next(iter(self.profiles.values())).plot(significance_level=alpha)
+            next(iter(self.profiles.values())).plot(
+                significance_level=alpha,
+                std_err=std_errs[0] if with_stderr else None,
+            )
         else:
             plot_idx = np.reshape(np.arange(1, len(self.profiles) + 1), (-1, 2)).T.flatten()
-            for idx, profile in zip(plot_idx, self.profiles.values()):
+            for par_idx, (idx, profile) in enumerate(zip(plot_idx, self.profiles.values())):
                 plt.subplot(self.n_components, 2, idx)
-                profile.plot(significance_level=alpha)
+                profile.plot(
+                    significance_level=alpha,
+                    std_err=std_errs[par_idx] if with_stderr else None,
+                    **kwargs,
+                )
 
 
 @dataclass(frozen=True)