diff --git a/README.md b/README.md
index a301c27..0ffec25 100644
--- a/README.md
+++ b/README.md
@@ -99,7 +99,7 @@ Additional dependencies can be installed using
     # To add scipy-related dependencies
     pip install continuous-timeseries[scipy]
 
-    # To add all optionatl dependencies
+    # To add all optional dependencies
     pip install continuous-timeseries[full]
     ```
 
diff --git a/changelog/18.docs.md b/changelog/18.docs.md
new file mode 100644
index 0000000..4f85fb1
--- /dev/null
+++ b/changelog/18.docs.md
@@ -0,0 +1 @@
+Added a tutorial into our support for creating emissions pathways that are compatible with a given budget (see [Budget-compatible emissions pathways](../tutorials/budget_compatible_pathways)).
diff --git a/changelog/18.feature.md b/changelog/18.feature.md
new file mode 100644
index 0000000..ec7cec3
--- /dev/null
+++ b/changelog/18.feature.md
@@ -0,0 +1 @@
+Added [`budget_compatible_pathways`][continuous_timeseries.budget_compatible_pathways] to support the creation of pathways compatible with a given budget.
diff --git a/changelog/18.trivial.md b/changelog/18.trivial.md
new file mode 100644
index 0000000..c37861c
--- /dev/null
+++ b/changelog/18.trivial.md
@@ -0,0 +1 @@
+Added support for rending maths with latex in the docs.
diff --git a/docs/NAVIGATION.md b/docs/NAVIGATION.md
index a663dea..6e9e638 100644
--- a/docs/NAVIGATION.md
+++ b/docs/NAVIGATION.md
@@ -10,6 +10,7 @@ See https://oprypin.github.io/mkdocs-literate-nav/
 - [Tutorials](tutorials/index.md)
     - [Timeseries tutorial](tutorials/timeseries_tutorial.py)
     - [Higher-order interpolation](tutorials/higher_order_interpolation.py)
+    - [Budget-compatible emissions pathways](tutorials/budget_compatible_pathways.py)
     - [Continuous timeseries tutorial](tutorials/continuous_timeseries_tutorial.py)
     - [Discrete timeseries tutorial](tutorials/discrete_timeseries_tutorial.py)
 - [Further background](further-background/index.md)
diff --git a/docs/javascripts/katex.js b/docs/javascripts/katex.js
new file mode 100644
index 0000000..f7fd704
--- /dev/null
+++ b/docs/javascripts/katex.js
@@ -0,0 +1,10 @@
+document$.subscribe(({ body }) => {
+  renderMathInElement(body, {
+    delimiters: [
+      { left: "$$",  right: "$$",  display: true },
+      { left: "$",   right: "$",   display: false },
+      { left: "\\(", right: "\\)", display: false },
+      { left: "\\[", right: "\\]", display: true }
+    ],
+  })
+})
diff --git a/docs/tutorials/budget_compatible_pathways.py b/docs/tutorials/budget_compatible_pathways.py
new file mode 100644
index 0000000..1e660f4
--- /dev/null
+++ b/docs/tutorials/budget_compatible_pathways.py
@@ -0,0 +1,396 @@
+# ---
+# jupyter:
+#   jupytext:
+#     text_representation:
+#       extension: .py
+#       format_name: percent
+#       format_version: '1.3'
+#       jupytext_version: 1.16.6
+#   kernelspec:
+#     display_name: Python 3 (ipykernel)
+#     language: python
+#     name: python3
+# ---
+
+# %% [markdown]
+# # Budget-compatible pathways
+#
+# Here we discuss our module
+# for creating emissions pathways compatible with a given budget.
+
+# %% [markdown]
+# ## Imports
+
+# %%
+import matplotlib.pyplot as plt
+import numpy as np
+import openscm_units
+import pint
+import scipy.interpolate
+
+import continuous_timeseries as ct
+import continuous_timeseries.budget_compatible_pathways as ct_bcp
+from continuous_timeseries.timeseries_continuous import ContinuousFunctionScipyPPoly
+
+# %% [markdown]
+# ## Set up pint
+
+# %%
+pint.set_application_registry(openscm_units.unit_registry)
+
+# %% [markdown]
+# ## Handy pint aliases
+
+# %%
+UR = pint.get_application_registry()
+Q = UR.Quantity
+
+# %% [markdown]
+# ## Set up matplotlib to work with pint
+#
+# For details, see the pint docs
+# ([stable docs](https://pint.readthedocs.io/en/stable/user/plotting.html),
+# [last version that we checked at the time of writing](https://pint.readthedocs.io/en/0.24.4/user/plotting.html))
+# [or our docs on unit-aware plotting](../discrete_timeseries_tutorial#unit-aware-plotting).
+
+# %%
+UR.setup_matplotlib(enable=True)
+
+# %% [markdown]
+# ## Calculating a budget-compatible pathway
+
+# %% [markdown]
+# Let's imagine we start with some emissions budget from some point in time.
+
+# %%
+budget = Q(500, "GtCO2")
+budget_start_time = Q(2020.0, "yr")
+
+# %% [markdown]
+# ### Linear pathway to zero
+#
+# Calculating a linear pathway to zero emissions
+# that is compatible with this budget is trivial.
+# The only other piece of information you need
+# is the emissions you want to start from.
+# Normally, that will be emissions at the time
+# from which the budget is/was available.
+
+# %%
+emissions_start = Q(10.4, "GtC / yr").to("GtCO2 / yr")
+
+# %%
+linear_in_line_with_budget = ct_bcp.derive_linear_path(
+    budget=budget,
+    budget_start_time=budget_start_time,
+    emissions_start=emissions_start,
+)
+linear_in_line_with_budget
+
+# %%
+fig, axes = plt.subplots(nrows=2, sharex=True)
+
+linear_in_line_with_budget.plot(ax=axes[0])
+linear_in_line_with_budget.integrate(Q(0, "GtC"), name_res="Cumulative emissions").plot(
+    ax=axes[1]
+)
+
+for ax in axes:
+    ax.legend()
+    ax.grid()
+
+fig.tight_layout()
+
+# %% [markdown]
+# There is one thing to notice on this plot.
+# The pathway is extended slightly beyond the net-zero year with zero emissions.
+# This ensures that there is at least one full year with zero emissions
+# and you get sensible results if you extrapolate after creating the pathway.
+
+# %%
+extrapolation_times = np.arange(
+    np.ceil(budget_start_time).to("yr").m, 2100.0 + 1.0, 1.0
+) * UR.Unit("yr")
+
+linear_extrapolated = linear_in_line_with_budget.interpolate(
+    extrapolation_times, allow_extrapolation=True
+)
+linear_extrapolated.timeseries_continuous.name = "Extrapolated"
+
+fig, axes = plt.subplots(nrows=2, sharex=True)
+
+
+linear_extrapolated.plot(ax=axes[0])
+linear_extrapolated.integrate(
+    Q(0, "GtC"), name_res="Extrapolated cumulative emissions"
+).plot(ax=axes[1])
+
+for ax in axes:
+    ax.legend()
+    ax.grid()
+
+fig.tight_layout()
+
+# %% [markdown]
+# #### Annual steps
+#
+# In general, countries don't report emissions instantaneously.
+# Instead, they report emissions as the average over each year.
+# With Continuous Timeseries, it is trivial to transform our pathway
+# into something which can be directly compared with country emissions.
+
+# %%
+linear_annual_steps_in_line_with_budget = ct_bcp.convert_to_annual_constant_emissions(
+    linear_in_line_with_budget, name_res="Annualised"
+)
+
+# %%
+continuous_plot_kwargs = dict(linewidth=2.0, alpha=0.7)
+pathways = [linear_in_line_with_budget, linear_annual_steps_in_line_with_budget]
+
+fig, axes = plt.subplots(nrows=2, sharex=True, figsize=(10, 5))
+
+for pw in pathways:
+    pw.plot(ax=axes[0], continuous_plot_kwargs=continuous_plot_kwargs)
+    pw.integrate(Q(0, "GtC"), name_res=f"Cumulative emissions {pw.name}").plot(
+        ax=axes[1], continuous_plot_kwargs=continuous_plot_kwargs
+    )
+
+for ax in axes:
+    ax.legend()
+    ax.grid()
+
+fig.tight_layout()
+
+# %% [markdown]
+# The difference between the two clearer if we take a more extreme example.
+
+# %%
+budget_small = Q(0.425, "GtCO2")
+budget_small_start_time = Q(2020.0, "yr")
+emissions_small_start = Q(300, "MtCO2 / yr")
+
+# %%
+linear_in_line_with_budget_small = ct_bcp.derive_linear_path(
+    budget=budget_small,
+    budget_start_time=budget_small_start_time,
+    emissions_start=emissions_small_start,
+)
+linear_annual_steps_in_line_with_budget_small = (
+    ct_bcp.convert_to_annual_constant_emissions(
+        linear_in_line_with_budget_small, name_res="Annualised"
+    )
+)
+
+# %%
+continuous_plot_kwargs = dict(linewidth=2.0, alpha=0.7)
+pathways = [
+    linear_in_line_with_budget_small,
+    linear_annual_steps_in_line_with_budget_small,
+]
+
+fig, axes = plt.subplots(nrows=2, sharex=True, figsize=(10, 5))
+
+for pw in pathways:
+    pw.plot(ax=axes[0], continuous_plot_kwargs=continuous_plot_kwargs)
+    pw.integrate(Q(0, "GtC"), name_res=f"Cumulative emissions {pw.name}").plot(
+        ax=axes[1], continuous_plot_kwargs=continuous_plot_kwargs
+    )
+
+for ax in axes:
+    ax.legend()
+    ax.grid()
+
+fig.tight_layout()
+
+# %% [markdown]
+# This example also shows why annualising emissions can be tricky.
+# The drop from one year to the next is the same for all steps,
+# except the year in which (instantaneous) net zero is reached,
+# which has a smaller drop to the zero emissions.
+# For example, in the annualised pathway above,
+# the drop in emissions from 2021 to 2022 is around 100 MtCO<sub>2</sub> / yr
+# while the drop in emissions from 2022 to 2023 is only around 40 MtCO<sub>2</sub> / yr.
+# (There are ways to do such calculations in their discrete form,
+# Continuous Timeseries just removes those headaches.)
+
+# %% [markdown]
+# ### Curved pathway to zero
+#
+# A linear pathway is one way to get to zero.
+# Other options are possible, for example a smoother pathway.
+#
+# #### Quadratic
+#
+# First we show a symmetric quadratic.
+# This generally produces nice results.
+# Its major downside is that its gradient is zero today,
+# which may not be desirable in all cases.
+
+# %%
+quadratic_in_line_with_budget = ct_bcp.derive_symmetric_quadratic_path(
+    budget=budget,
+    budget_start_time=budget_start_time,
+    emissions_start=emissions_start,
+)
+
+# %%
+continuous_plot_kwargs = dict(linewidth=2.0, alpha=0.7)
+pathways = [linear_in_line_with_budget, quadratic_in_line_with_budget]
+
+fig, axes = plt.subplots(nrows=2, sharex=True, figsize=(10, 5))
+
+for pw in pathways:
+    pw.plot(ax=axes[0], continuous_plot_kwargs=continuous_plot_kwargs)
+    pw.integrate(Q(0, "GtC"), name_res=f"Cumulative emissions {pw.name}").plot(
+        ax=axes[1], continuous_plot_kwargs=continuous_plot_kwargs
+    )
+
+for ax in axes:
+    ax.legend()
+    ax.grid()
+
+fig.tight_layout()
+
+# %%
+quadratic_annual_steps_in_line_with_budget = (
+    ct_bcp.convert_to_annual_constant_emissions(
+        quadratic_in_line_with_budget, name_res="Annualised quadratic"
+    )
+)
+
+# %%
+continuous_plot_kwargs = dict(linewidth=2.0, alpha=0.7)
+pathways = [quadratic_in_line_with_budget, quadratic_annual_steps_in_line_with_budget]
+
+fig, axes = plt.subplots(nrows=2, sharex=True, figsize=(10, 5))
+
+for pw in pathways:
+    pw.plot(ax=axes[0], continuous_plot_kwargs=continuous_plot_kwargs)
+    pw.integrate(Q(0, "GtC"), name_res=f"Cumulative emissions {pw.name}").plot(
+        ax=axes[1], continuous_plot_kwargs=continuous_plot_kwargs
+    )
+
+for ax in axes:
+    ax.legend()
+    ax.grid()
+
+fig.tight_layout()
+
+
+# %% [markdown]
+# #### Custom
+#
+# Here we show how to make a custom fit,
+# in case you want to explore further yourself.
+# Here we do the fit such that the emissions and gradient today
+# are specified, but we still stay within the budget.
+# We do this by fitting a cubic pathway.
+
+
+# %%
+def get_cubic_pathway(  # noqa: D103, PLR0913
+    budget,
+    budget_start_time,
+    emissions_start,
+    emissions_gradient_start,
+    name_res,
+    max_iter: int = 100,
+):
+    time_units = budget_start_time.u
+    values_units = emissions_start.u
+
+    E_0 = emissions_start
+    m_0 = emissions_gradient_start
+
+    # A first guess, doesn't really matter what it is.
+    # (Would probably break in some cases, but not the concern here).
+    nz_yr = budget_start_time + Q(1.0, "yr")
+    time_bounds = np.hstack([budget_start_time, nz_yr])
+    # Figure out scale factor to match the budget
+    for _ in range(max_iter):
+        # Handy constant
+        nzd = nz_yr - budget_start_time
+
+        b = -(2 * m_0 * nzd + 3 * E_0) / (nzd**2)
+        a = (m_0 * nzd + 2 * E_0) / (nzd**3)
+
+        c_non_zero = np.array(
+            [
+                [a.to(values_units / time_units**3).m],
+                [b.to(values_units / time_units**2).m],
+                [m_0.to(values_units / time_units).m],
+                [E_0.to(values_units).m],
+            ]
+        )
+
+        ppoly_raw = scipy.interpolate.PPoly(x=time_bounds.m, c=c_non_zero)
+        budget_raw = ppoly_raw.integrate(
+            budget_start_time.to(time_units).m, nz_yr.to(time_units).m
+        )
+        scale_factor = budget.to(values_units * time_units).m / budget_raw
+        if np.isclose(scale_factor, 1.0, rtol=1e-4):
+            break
+
+        nz_yr = budget_start_time + nzd * scale_factor
+        time_bounds = np.hstack([budget_start_time, nz_yr])
+
+    else:
+        msg = "Did not converge"
+        raise AssertionError(msg)
+
+    # Add on the zero component
+    c = np.hstack([c_non_zero, np.zeros((c_non_zero.shape[0], 1))])
+    last_ts_time = np.floor(nz_yr) + 2.0 * nz_yr.to("yr").u
+    x_bounds = np.hstack([budget_start_time, nz_yr, last_ts_time])
+    x = x_bounds.to(time_units).m
+
+    ppoly = scipy.interpolate.PPoly(x=x, c=c)
+
+    tsc = ct.TimeseriesContinuous(
+        name=name_res,
+        time_units=time_units,
+        values_units=values_units,
+        function=ContinuousFunctionScipyPPoly(ppoly),
+        domain=(budget_start_time, last_ts_time),
+    )
+    ts = ct.Timeseries(
+        time_axis=ct.TimeAxis(x_bounds),
+        timeseries_continuous=tsc,
+    )
+
+    return ts
+
+
+# %%
+emissions_gradient_start = Q(-1.0, "GtCO2 / yr / yr")
+times_plot = (
+    np.arange(budget_start_time.to("yr").m, 2081, 1.0) * budget_start_time.to("yr").u
+)
+
+fig, axes = plt.subplots(nrows=3, sharex=True, figsize=(10, 6))
+for emissions_gradient_start in Q([-2.0, 0.0, 5.0], "GtCO2 / yr / yr"):
+    cubic = get_cubic_pathway(
+        budget=budget,
+        budget_start_time=budget_start_time,
+        emissions_start=emissions_start,
+        emissions_gradient_start=emissions_gradient_start,
+        name_res=f"Initial gradient: {emissions_gradient_start:.2f}",
+    ).interpolate(times_plot, allow_extrapolation=True)
+
+    annualised = ct_bcp.convert_to_annual_constant_emissions(
+        cubic, name_res=f"Annualised {cubic.name}"
+    )
+
+    cubic.plot(ax=axes[0], continuous_plot_kwargs=dict(linewidth=2.0))
+    annualised.plot(
+        ax=axes[0], continuous_plot_kwargs=dict(label="", zorder=1.0, alpha=0.6)
+    )
+    cubic.differentiate().plot(ax=axes[1])
+    cubic.integrate(Q(0, "GtC")).plot(ax=axes[2])
+
+axes[0].legend()
+for ax in axes:
+    ax.grid()
+
+fig.tight_layout()
diff --git a/docs/tutorials/index.md b/docs/tutorials/index.md
index a96f22c..e2386bd 100644
--- a/docs/tutorials/index.md
+++ b/docs/tutorials/index.md
@@ -29,6 +29,12 @@ See [higher-order interpolation](./higher_order_interpolation)
 for a tutorial on how to do higher-order interpolation with Continuous Timeseries
 and some tips on how to go beyond the default options available.
 
+## Budget-compatible emissions pathways
+
+See [budget-compatible emissions pathways](./budget_compatible_pathways)
+for a tutorial on how to calculate emissions pathways
+compatible with a given budget.
+
 ## Continuous timeseries handling
 
 See [Continuous timeseries](./continuous_timeseries_tutorial)
diff --git a/mkdocs.yml b/mkdocs.yml
index f9973ff..eaa590a 100644
--- a/mkdocs.yml
+++ b/mkdocs.yml
@@ -110,6 +110,9 @@ plugins:
 markdown_extensions:
   # https://squidfunk.github.io/mkdocs-material/setup/extensions/python-markdown/#attribute-lists
   - attr_list
+  # https://squidfunk.github.io/mkdocs-material/reference/math/#katex-mkdocsyml
+  - pymdownx.arithmatex:
+      generic: true
   # Allow admonitions, useful for deprecation warnings
   # https://facelessuser.github.io/pymdown-extensions/extensions/blocks/plugins/admonition/
   - pymdownx.blocks.admonition
@@ -138,6 +141,14 @@ markdown_extensions:
   - toc:
       permalink: "#"
 
+extra_javascript:
+  - javascripts/katex.js
+  - https://unpkg.com/katex@0/dist/katex.min.js
+  - https://unpkg.com/katex@0/dist/contrib/auto-render.min.js
+
+extra_css:
+  - https://unpkg.com/katex@0/dist/katex.min.css
+
 watch:
   - README.md
   # Auto-generate if `src` changes (because this changes API docs)
diff --git a/pdm.lock b/pdm.lock
index 3f2105a..954d374 100644
--- a/pdm.lock
+++ b/pdm.lock
@@ -5,14 +5,14 @@
 groups = ["default", "all-dev", "dev", "docs", "full", "plots", "scipy", "tests"]
 strategy = ["inherit_metadata"]
 lock_version = "4.5.0"
-content_hash = "sha256:68f13428514726753e6dd3e970ee62c0dfd4cbc731f52a04e93e2f1e1fd4306f"
+content_hash = "sha256:d588a6199c537a80be2873d5863598ca348c864dc25da7316cb96892c8df50dc"
 
 [[metadata.targets]]
 requires_python = ">=3.9"
 
 [[package]]
 name = "anyio"
-version = "4.7.0"
+version = "4.8.0"
 requires_python = ">=3.9"
 summary = "High level compatibility layer for multiple asynchronous event loop implementations"
 groups = ["all-dev", "docs"]
@@ -23,8 +23,8 @@ dependencies = [
     "typing-extensions>=4.5; python_version < \"3.13\"",
 ]
 files = [
-    {file = "anyio-4.7.0-py3-none-any.whl", hash = "sha256:ea60c3723ab42ba6fff7e8ccb0488c898ec538ff4df1f1d5e642c3601d07e352"},
-    {file = "anyio-4.7.0.tar.gz", hash = "sha256:2f834749c602966b7d456a7567cafcb309f96482b5081d14ac93ccd457f9dd48"},
+    {file = "anyio-4.8.0-py3-none-any.whl", hash = "sha256:b5011f270ab5eb0abf13385f851315585cc37ef330dd88e27ec3d34d651fd47a"},
+    {file = "anyio-4.8.0.tar.gz", hash = "sha256:1d9fe889df5212298c0c0723fa20479d1b94883a2df44bd3897aa91083316f7a"},
 ]
 
 [[package]]
@@ -726,7 +726,7 @@ name = "flexcache"
 version = "0.3"
 requires_python = ">=3.9"
 summary = "Saves and loads to the cache a transformed versions of a source object."
-groups = ["default", "all-dev", "docs"]
+groups = ["default", "all-dev", "docs", "tests"]
 dependencies = [
     "typing-extensions",
 ]
@@ -740,7 +740,7 @@ name = "flexparser"
 version = "0.4"
 requires_python = ">=3.9"
 summary = "Parsing made fun ... using typing."
-groups = ["default", "all-dev", "docs"]
+groups = ["default", "all-dev", "docs", "tests"]
 dependencies = [
     "typing-extensions",
 ]
@@ -831,7 +831,7 @@ files = [
 name = "globalwarmingpotentials"
 version = "0.11.1"
 summary = "Global warming potentials of greenhouse gases from various IPCC reports"
-groups = ["all-dev", "docs"]
+groups = ["all-dev", "docs", "tests"]
 files = [
     {file = "globalwarmingpotentials-0.11.1-py2.py3-none-any.whl", hash = "sha256:0f208afdaea1e464b1d2f760f06ab09e3870876c48a322252536809a5087ed11"},
     {file = "globalwarmingpotentials-0.11.1.tar.gz", hash = "sha256:c0d98d6fe7fb7d21bc2c17d7b41dfdf5a08b3a80778943e1f0d2f7f29372038b"},
@@ -2108,7 +2108,7 @@ name = "openscm-units"
 version = "0.6.3"
 requires_python = ">=3.9"
 summary = "Handling of units related to simple climate modelling."
-groups = ["all-dev", "docs"]
+groups = ["all-dev", "docs", "tests"]
 dependencies = [
     "globalwarmingpotentials",
     "numpy",
@@ -2175,7 +2175,7 @@ name = "pandas"
 version = "2.2.3"
 requires_python = ">=3.9"
 summary = "Powerful data structures for data analysis, time series, and statistics"
-groups = ["all-dev", "docs"]
+groups = ["all-dev", "docs", "tests"]
 dependencies = [
     "numpy>=1.22.4; python_version < \"3.11\"",
     "numpy>=1.23.2; python_version == \"3.11\"",
@@ -2361,7 +2361,7 @@ name = "pint"
 version = "0.24.4"
 requires_python = ">=3.9"
 summary = "Physical quantities module"
-groups = ["default", "all-dev", "docs"]
+groups = ["default", "all-dev", "docs", "tests"]
 dependencies = [
     "flexcache>=0.3",
     "flexparser>=0.4",
@@ -2389,7 +2389,7 @@ name = "platformdirs"
 version = "4.3.6"
 requires_python = ">=3.8"
 summary = "A small Python package for determining appropriate platform-specific dirs, e.g. a `user data dir`."
-groups = ["default", "all-dev", "dev", "docs"]
+groups = ["default", "all-dev", "dev", "docs", "tests"]
 files = [
     {file = "platformdirs-4.3.6-py3-none-any.whl", hash = "sha256:73e575e1408ab8103900836b97580d5307456908a03e92031bab39e4554cc3fb"},
     {file = "platformdirs-4.3.6.tar.gz", hash = "sha256:357fb2acbc885b0419afd3ce3ed34564c13c9b95c89360cd9563f73aa5e2b907"},
@@ -2500,13 +2500,13 @@ files = [
 
 [[package]]
 name = "pygments"
-version = "2.18.0"
+version = "2.19.0"
 requires_python = ">=3.8"
 summary = "Pygments is a syntax highlighting package written in Python."
 groups = ["all-dev", "dev", "docs", "tests"]
 files = [
-    {file = "pygments-2.18.0-py3-none-any.whl", hash = "sha256:b8e6aca0523f3ab76fee51799c488e38782ac06eafcf95e7ba832985c8e7b13a"},
-    {file = "pygments-2.18.0.tar.gz", hash = "sha256:786ff802f32e91311bff3889f6e9a86e81505fe99f2735bb6d60ae0c5004f199"},
+    {file = "pygments-2.19.0-py3-none-any.whl", hash = "sha256:4755e6e64d22161d5b61432c0600c923c5927214e7c956e31c23923c89251a9b"},
+    {file = "pygments-2.19.0.tar.gz", hash = "sha256:afc4146269910d4bdfabcd27c24923137a74d562a23a320a41a55ad303e19783"},
 ]
 
 [[package]]
@@ -2604,7 +2604,7 @@ name = "python-dateutil"
 version = "2.9.0.post0"
 requires_python = "!=3.0.*,!=3.1.*,!=3.2.*,>=2.7"
 summary = "Extensions to the standard Python datetime module"
-groups = ["all-dev", "docs", "full", "plots"]
+groups = ["all-dev", "docs", "full", "plots", "tests"]
 dependencies = [
     "six>=1.5",
 ]
@@ -2631,7 +2631,7 @@ files = [
 name = "pytz"
 version = "2024.2"
 summary = "World timezone definitions, modern and historical"
-groups = ["all-dev", "docs"]
+groups = ["all-dev", "docs", "tests"]
 files = [
     {file = "pytz-2024.2-py2.py3-none-any.whl", hash = "sha256:31c7c1817eb7fae7ca4b8c7ee50c72f93aa2dd863de768e1ef4245d426aa0725"},
     {file = "pytz-2024.2.tar.gz", hash = "sha256:2aa355083c50a0f93fa581709deac0c9ad65cca8a9e9beac660adcbd493c798a"},
@@ -3210,7 +3210,7 @@ name = "six"
 version = "1.17.0"
 requires_python = "!=3.0.*,!=3.1.*,!=3.2.*,>=2.7"
 summary = "Python 2 and 3 compatibility utilities"
-groups = ["all-dev", "docs", "full", "plots"]
+groups = ["all-dev", "docs", "full", "plots", "tests"]
 files = [
     {file = "six-1.17.0-py2.py3-none-any.whl", hash = "sha256:4721f391ed90541fddacab5acf947aa0d3dc7d27b2e1e8eda2be8970586c3274"},
     {file = "six-1.17.0.tar.gz", hash = "sha256:ff70335d468e7eb6ec65b95b99d3a2836546063f63acc5171de367e834932a81"},
@@ -3412,7 +3412,7 @@ name = "tzdata"
 version = "2024.2"
 requires_python = ">=2"
 summary = "Provider of IANA time zone data"
-groups = ["all-dev", "docs"]
+groups = ["all-dev", "docs", "tests"]
 files = [
     {file = "tzdata-2024.2-py2.py3-none-any.whl", hash = "sha256:a48093786cdcde33cad18c2555e8532f34422074448fbc874186f0abd79565cd"},
     {file = "tzdata-2024.2.tar.gz", hash = "sha256:7d85cc416e9382e69095b7bdf4afd9e3880418a2413feec7069d533d6b4e31cc"},
diff --git a/pyproject.toml b/pyproject.toml
index 897d9fd..065546f 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -68,6 +68,7 @@ docs = [
 tests = [
     "pytest>=8.3.3",
     "coverage>=7.6.0",
+    "openscm-units>=0.6.3",
     "pytest-cov>=5.0.0",
     "pytest-regressions>=2.6.0",
     "ipython>=8.18.1",
@@ -185,6 +186,9 @@ ignore = [
 "docs/how-to-guides/how-to-make-a-step-forcing.py" = [
     "E501",
 ]
+"docs/tutorials/budget_compatible_pathways.py" = [
+    "E501",
+]
 "docs/tutorials/continuous_timeseries_tutorial.py" = [
     "E501",
 ]
diff --git a/src/continuous_timeseries/budget_compatible_pathways.py b/src/continuous_timeseries/budget_compatible_pathways.py
new file mode 100644
index 0000000..76d6578
--- /dev/null
+++ b/src/continuous_timeseries/budget_compatible_pathways.py
@@ -0,0 +1,501 @@
+"""
+Creation of emissions pathways compatible with a given budget
+"""
+
+from __future__ import annotations
+
+import numpy as np
+
+from continuous_timeseries.discrete_to_continuous import InterpolationOption
+from continuous_timeseries.exceptions import MissingOptionalDependencyError
+from continuous_timeseries.time_axis import TimeAxis
+from continuous_timeseries.timeseries import Timeseries
+from continuous_timeseries.timeseries_continuous import (
+    ContinuousFunctionScipyPPoly,
+    TimeseriesContinuous,
+)
+from continuous_timeseries.typing import PINT_NUMPY_ARRAY, PINT_SCALAR
+
+
+def calculate_linear_net_zero_time(
+    budget: PINT_SCALAR,
+    budget_start_time: PINT_SCALAR,
+    emissions_start: PINT_SCALAR,
+) -> PINT_SCALAR:
+    """
+    Calculate the net-zero time, assuming a linear path to net-zero
+
+    Parameters
+    ----------
+    budget
+        Budget to match
+
+    budget_start_time
+        Time from which the budget is available.
+
+        E.g., if the budget is available from 1 Jan 2025,
+        supply `Q(2025.0, "yr")` (assuming you're working with pint quantities).
+
+    emissions_start
+        Emissions from which the linear path should start
+
+    Returns
+    -------
+    :
+        Net-zero time
+    """
+    net_zero_time: PINT_SCALAR = budget_start_time + 2 * budget / emissions_start
+
+    return net_zero_time
+
+
+def derive_linear_path(
+    budget: PINT_SCALAR,
+    budget_start_time: PINT_SCALAR,
+    emissions_start: PINT_SCALAR,
+    name_res: str | None = None,
+) -> Timeseries:
+    r"""
+    Derive a linear pathway that stays within a given budget
+
+    Parameters
+    ----------
+    budget
+        Budget to match
+
+    budget_start_time
+        Time from which the budget is available.
+
+        E.g., if the budget is available from 1 Jan 2025,
+        supply `Q(2025.0, "yr")` (assuming you're working with pint quantities).
+
+    emissions_start
+        Emissions from which the linear path should start
+
+    name_res
+        Name to use for the result.
+
+        If not supplied, we use a verbose but clear default.
+
+    Returns
+    -------
+    :
+        Linear pathway to net-zero in line with the budget
+
+    Notes
+    -----
+    We're solving for emissions, $y$, as a function of time, $x$.
+    We pick a linear emissions pathway
+
+    $$
+    y(x) = e_0 * (1 - \frac{x - x_0}{x_{nz} - x_0})
+    $$
+
+    Simplifying slightly, we have
+
+    $$
+    y(x) = e_0 * \frac{x_{nz} - x}{x_{nz} - x_0}
+    $$
+
+    where $e_0$ is emissions at the known time (normally today), $x_0$,
+    and $x_{nz}$ is the net-zero time.
+
+    By geometry, the integral of this curve between $x_0$ and $x_nz$ is:
+
+    $$
+    \frac{1}{2} (x_0 - x_{nz}) * e_0
+    $$
+
+    You can also do this with calculus:
+
+    $$
+    \begin{equation}
+    \begin{split}
+    \int_{x_0}^{x_{nz}} y(x) \, dx &= \int_{x_0}^{x_{nz}} e_0 * (x_{nz} - x) / (x_{nz} - x_0) \, dx \\
+                             &= [-e_0 (x_{nz} - x)^2 / (2 * (x_{nz} - x_0))]_{x_0}^{x_{nz}} \\
+                             &= -e_0 (x_{nz} - x_0)^2 / (2 * (x_{nz} - x_0))) - e_0 (x_{nz} - x_{nz})^2 / (2 * (x_{nz} - x_0))) \\
+                             &= e_0 (x_0 - x_{nz}) / 2
+    \end{split}
+    \end{equation}
+    $$
+
+    This integral should be equal to the allowed buget:
+
+    $$
+    \frac{1}{2} (x_0 - x_{nz}) * e_0 = \text{budget}
+    $$
+
+    therefore
+
+    $$
+    x_{nz} = x_0 + \frac{2 * \text{budget}}{e_0}
+    $$
+    """  # noqa: E501
+    if name_res is None:
+        name_res = (
+            "Linear emissions\n"
+            f"compatible with budget of {budget:.2f}\n"
+            f"from {budget_start_time:.2f}"
+        )
+
+    net_zero_time = calculate_linear_net_zero_time(
+        budget=budget,
+        budget_start_time=budget_start_time,
+        emissions_start=emissions_start,
+    )
+
+    last_ts_time = np.floor(net_zero_time) + 2.0 * net_zero_time.to("yr").u
+
+    time_axis_bounds: PINT_NUMPY_ARRAY = np.hstack(  # type: ignore # mypy confused by pint
+        [budget_start_time, net_zero_time, last_ts_time]
+    )
+    values_at_bounds: PINT_NUMPY_ARRAY = np.hstack(  # type: ignore # mypy confused by pint
+        [emissions_start, 0.0 * emissions_start, 0.0 * emissions_start]
+    )
+
+    emms_linear_pathway = Timeseries.from_arrays(
+        time_axis_bounds=time_axis_bounds,
+        values_at_bounds=values_at_bounds,
+        interpolation=InterpolationOption.Linear,
+        name=name_res,
+    )
+
+    return emms_linear_pathway
+
+
+def derive_symmetric_quadratic_path(
+    budget: PINT_SCALAR,
+    budget_start_time: PINT_SCALAR,
+    emissions_start: PINT_SCALAR,
+    name_res: str | None = None,
+) -> Timeseries:
+    r"""
+    Derive a quadratic pathway that stays within a given budget
+
+    The major downside of this approach is
+    that the gradient of the output path is zero at `budget_start_time`.
+
+    Parameters
+    ----------
+    budget
+        Budget to match
+
+    budget_start_time
+        Time from which the budget is available.
+
+        E.g., if the budget is available from 1 Jan 2025,
+        supply `Q(2025.0, "yr")` (assuming you're working with pint quantities).
+
+    emissions_start
+        Emissions from which the quadratic path should start.
+
+    name_res
+        Name to use for the result.
+
+        If not supplied, we use a verbose but clear default.
+
+    Returns
+    -------
+    :
+        Symmetric quadratic pathway to net-zero in line with the budget
+
+    Notes
+    -----
+    We're solving for emissions, $y$, as a function of time, $x$.
+    We pick a quadratic emissions pathway comprised of two pieces.
+    The two pieces are equal in size and split the time period
+    from the known time (normally today), $x_0$,
+    to the net zero time, $x_{nz}$, in two.
+    For convenience, we also define the halfway to net-zero point,
+    $x_{nzh} = \frac{x_0 + x_{nz}}{2}$.
+
+    $$
+    y(x) = \begin{cases}
+       a_1 (x - x_0)^2 + b_1 (x - x_0) + c_1 &\text{if } x \leq x_{nzh} \\
+       a_2 (x - x_{nzh})^2 + b_2 (x - x_{nzh}) + c_2 &\text{otherwise}
+    \end{cases}
+    $$
+
+    Therefore
+
+    $$
+    \frac{dy}{dx}(x) = \begin{cases}
+       2 a_1 (x - x_0) + b_1 &\text{if } x \leq x_{nzh} \\
+       2 a_2 (x - x_{nzh}) + b_2 &\text{otherwise}
+    \end{cases}
+    $$
+
+    To set the constants, we need some boundary conditions.
+
+    At $x = x_0$, emissions should be equal to the starting emissions, $e_0$.
+    This immediately sets $c_1 = e_0$.
+
+    We also choose to set the quadratics such that
+    at the time point which is halfway to net zero,
+    emissions are half of their original value.
+    Thus, at $x = x_{nzh}$, emissions should be $e_0 / 2$.
+    This immediately sets $c_2 = e_0 / 2$.
+
+    This condition also means that the decrease in emissions should be the same
+    between $x_0$ and $x_{nzh}$ as it is between $x_{nzh}$ and $x_{nz}$.
+    We also want the gradient to be continuous
+    at the boundary between the two quadratics.
+
+    As a result, we can see that the two quadratics are simply
+    a translation and a reflection of each other
+    (they have the same change between their defining points
+    and the same gradient at the boundary between the two intervals,
+    so there are no more degrees of freedom with which we could
+    introduce different shapes),
+    i.e. they are symmetric (about a carefully chosen axis).
+
+    By the symmetry argument, we have that $a_1 = -a_2$
+    and the gradient at $x=x_0$ must equal the gradient at $x=x_{nx}$.
+    The gradient at $x=x_{nx}$ is zero, therefore
+
+    $$
+    \frac{dy}{dx}(x_0) = 2 a_1 (x_0 - x_0) + b_1 = 0.0 \Rightarrow b_1 = 0.0
+    $$
+
+    We can now use the fact that $y(x_{nzh}) = e_0 / 2$ to solve for $a_1$,
+    and therefore also $a_2$.
+
+    $$
+    \begin{equation}
+    \begin{split}
+    y(x_{nzh}) &= e_0 / 2 \\
+    a_1 (x_{nzh} - x_0)^2 + e_0 &= e_0 / 2 \\
+    a_1 &= -\frac{e_0}{2 (x_{nzh} - x_0)^2}
+    \end{split}
+    \end{equation}
+    $$
+
+    The last constant to solve for is $b_2$.
+    We solve for this using the first-order continuity constraint at the boundary.
+
+    $$
+    \frac{dy}{dx}(x_{nzh}) = 2 a_1 (x_{nzh} - x_0) = 2 a_2 (x - x_{nzh}) + b_2 \Rightarrow b_2 = 2 a_1 (x_{nzh} - x_0)
+    $$
+
+    In summary, our constants are
+
+    $$
+    \begin{equation}
+    \begin{split}
+    a_1 &= -\frac{e_0}{2 (x_{nzh} - x_0)^2} \\
+    a_2 &= -a_1 \\
+    b_1 &= 0.0 \\
+    b_2 &= 2 a_1 (x_{nzh} - x_0) \\
+    c_1 &= e_0 \\
+    c_2 &= e_0 / 2
+    \end{split}
+    \end{equation}
+    $$
+
+    The last question is, where does the net-zero year come from?
+    To answer this, first consider the integral of our function
+    from $x_0$ to $x_{nz}$
+
+    $$
+    \begin{equation}
+    \begin{split}
+    \int_{x_0}^{x_{nz}} y(x) \, dx &= \int_{x_0}^{x_{nzh}} a_1 (x - x_0)^2 + c_1 \, dx \\
+                                    & \quad + \int_{x_{nzh}}^{x_{nz}} a_2 (x - x_{nzh})^2 + b_2 (x - x_{nzh}) + c_2 \, dx \\
+                                   &= [ a_1 (x - x_0)^3 / 3 + c_1 x ]_{x_0}^{x_{nzh}} \\
+                                    & + [ a_2 (x - x_{nzh})^3 / 3 + b_2 (x - x_{nzh})^2 / 2 + c_2 x ]_{x_{nzh}}^{x_{nz}} \\
+                                   &= a_1 (x_{nzh} - x_0)^3 / 3 + c_1 (x_{nzh} - x_0) \\
+                                    & + a_2 (x_{nz} - x_{nzh})^3 / 3 + b_2 (x_{nz} - x_{nzh})^2 / 2 + c_2 (x_{nz} - x_{nzh})
+    \end{split}
+    \end{equation}
+    $$
+
+    We next note that
+
+    $$
+    x_{nzh} - x_0 = \frac{x_0 + x_{nz}}{2} - x_0 = \frac{x_{nz} - x_0}{2} = \frac{2 x_{nz} - (x_{nz} + x_0)}{2} = x_{nz} - x_{nzh}
+    $$
+
+    Hence the cubic terms cancel because $a_1 = -a_2$ and we are left with
+
+    $$
+    \begin{equation}
+    \begin{split}
+    \int_{x_0}^{x_{nz}} y(x) \, dx &= c_1 (x_{nzh} - x_0) + b_2 (x_{nz} - x_{nzh})^2 / 2 + c_2 (x_{nzh} - x_0) \\
+                                   &= e_0 (x_{nzh} - x_0) + 2 a_1 (x_{nzh} - x_0) (x_{nz} - x_{nzh})^2 / 2 + e_0 / 2 (x_{nzh} - x_0) \\
+                                   &= \frac{3}{2} e_0 (x_{nzh} - x_0) - 2 \frac{e_0}{2 (x_{nzh} - x_0)^2} (x_{nzh} - x_0) (x_{nz} - x_{nzh})^2 / 2 \\
+                                   &= \frac{3}{2} e_0 (x_{nzh} - x_0) - \frac{e_0}{2 (x_{nzh} - x_0)} (x_{nz} - x_{nzh})^2 \\
+                                   &= \frac{3}{2} e_0 (x_{nzh} - x_0) - \frac{e_0}{2 (x_{nzh} - x_0)} (x_{nzh} - x_0)^2 \\
+                                   &= \frac{3}{2} e_0 (x_{nzh} - x_0) - \frac{e_0}{2} (x_{nzh} - x_0) \\
+                                   &= \frac{3}{2} e_0 (x_{nzh} - x_0) - \frac{e_0}{2} (x_{nzh} - x_0) \\
+                                   &= e_0 (x_{nzh} - x_0) \\
+                                   &= e_0 (\frac{x_{nz} - x_0}{2}) \\
+                                   &= \frac{1}{2} e_0 (x_{nz} - x_0)
+    \end{split}
+    \end{equation}
+    $$
+
+    Put more simply,
+    the integral is simply equal to the integral of a straight-line to net-zero.
+    Hence, we can simply use the net-zero year of a linear pathway to net-zero
+    in line with the budget,
+    and our quadratic pathway will have the same cumulative emissions
+    (i.e. will also match the budget).
+
+    As a result, our recipe is:
+
+    1. Calculate the net-zero year of a straight-line pathway to net-zero year.
+    1. Use that net-zero year to calculate our constants.
+    """  # noqa: E501
+    try:
+        import scipy.interpolate
+    except ImportError as exc:
+        raise MissingOptionalDependencyError(
+            "derive_symmetric_quadratic_path", requirement="scipy"
+        ) from exc
+
+    if name_res is None:
+        name_res = (
+            "Symmetric quadratic emissions\n"
+            f"compatible with budget of {budget:.2f}\n"
+            f"from {budget_start_time:.2f}"
+        )
+
+    time_units = budget_start_time.u
+    values_units = emissions_start.u
+
+    net_zero_time = calculate_linear_net_zero_time(
+        budget=budget,
+        budget_start_time=budget_start_time,
+        emissions_start=emissions_start,
+    )
+
+    e_0 = emissions_start
+    x_0 = budget_start_time
+    x_nz = net_zero_time
+    x_nzh = (x_0 + x_nz) / 2.0
+
+    a_1 = -e_0 / (2.0 * (x_nzh - x_0) ** 2)
+    a_2 = -a_1
+    b_1 = 0.0 * emissions_start.u / budget_start_time.u
+    b_2 = 2.0 * a_1 * (x_nzh - x_0)
+    c_1 = e_0
+    c_2 = e_0 / 2.0
+
+    a_coeffs: PINT_NUMPY_ARRAY = np.hstack([a_1, a_2])  # type: ignore # mypy confused ay pint
+    b_coeffs: PINT_NUMPY_ARRAY = np.hstack([b_1, b_2])  # type: ignore # mypy confused by pint
+    const_terms: PINT_NUMPY_ARRAY = np.hstack([c_1, c_2])  # type: ignore # mypy confused by pint
+
+    c_non_zero: PINT_NUMPY_ARRAY = np.vstack(  # type: ignore # mypy confused by pint
+        [
+            a_coeffs.to(values_units / time_units**2).m,
+            b_coeffs.to(values_units / time_units).m,
+            const_terms.to(values_units).m,
+        ]
+    )
+
+    # Ensure we have a zero tail
+    last_ts_time = np.floor(net_zero_time) + 2.0 * net_zero_time.to("yr").u
+    time_axis_bounds: PINT_NUMPY_ARRAY = np.hstack(  # type: ignore # mypy confused by pint
+        [x_0, x_nzh, x_nz, last_ts_time]
+    )
+
+    x = time_axis_bounds.to(net_zero_time.u).m
+    c = np.hstack([c_non_zero, np.zeros((c_non_zero.shape[0], 1))])
+
+    ppoly = scipy.interpolate.PPoly(c=c, x=x)
+    tsc = TimeseriesContinuous(
+        name=name_res,
+        time_units=time_units,
+        values_units=values_units,
+        function=ContinuousFunctionScipyPPoly(ppoly),
+        domain=(time_axis_bounds.min(), time_axis_bounds.max()),
+    )
+    emms_quadratic_pathway = Timeseries(
+        time_axis=TimeAxis(time_axis_bounds),
+        timeseries_continuous=tsc,
+    )
+
+    return emms_quadratic_pathway
+
+
+def convert_to_annual_time_axis(ts: Timeseries) -> Timeseries:
+    """
+    Convert a timeseries to an annual time axis
+
+    This is just a convenience method.
+    It has minimal checks, so may not always produce sensible results.
+
+    Parameters
+    ----------
+    ts
+        Timeseries to convert to an annual time axis.
+
+
+    Returns
+    -------
+    :
+        `ts`, with its time axis updated so that it has annual steps.
+
+        If `ts`'s first time point is not an integer,
+        it is also included in the output time axis
+        to ensure that `ts`'s cumulative emissions are conserved.
+    """
+    annual_time_axis = (
+        np.union1d(
+            ts.time_axis.bounds.min().to("yr").m,
+            np.arange(
+                np.ceil(ts.time_axis.bounds.min()).to("yr").m,
+                np.ceil(ts.time_axis.bounds.max()).to("yr").m + 1,
+                1.0,
+            ),
+        )
+        * ts.time_axis.bounds[0].to("yr").u
+    )
+
+    res = ts.interpolate(annual_time_axis, allow_extrapolation=True)
+
+    return res
+
+
+def convert_to_annual_constant_emissions(
+    ts: Timeseries, name_res: str | None = None
+) -> Timeseries:
+    """
+    Convert a timeseries to annual constant emissions
+
+    In other words, to annual-average emissions
+    (or annual-total, depending on how you think about emissions),
+    like what countries report.
+
+    If the time axis of `ts` starts with an integer year,
+    then you can simply sum the output emissions and you will get
+    the same cumulative emissions as `ts`.
+    If the time axis of `ts` does not start with an integer year,
+    then it is more complicated because your first time step
+    will not be a full year.
+
+    Parameters
+    ----------
+    ts
+        Timeseries to convert to an annual time axis.
+
+    name_res
+        Name to use for the result.
+
+        If not supplied, we use `f"{ts.name}_annualised"`.
+
+    Returns
+    -------
+    :
+        `ts`, with its time axis updated so that it has annual steps
+        and is piecewise constant.
+    """
+    if name_res is None:
+        name_res = f"{ts.name}_annualised"
+
+    annual_interp = convert_to_annual_time_axis(ts)
+    res = annual_interp.update_interpolation_integral_preserving(
+        InterpolationOption.PiecewiseConstantNextLeftClosed, name_res=name_res
+    )
+
+    return res
diff --git a/tests/integration/test_budget_compatible_pathways_integration.py b/tests/integration/test_budget_compatible_pathways_integration.py
new file mode 100644
index 0000000..c861d56
--- /dev/null
+++ b/tests/integration/test_budget_compatible_pathways_integration.py
@@ -0,0 +1,357 @@
+"""
+Integration tests of `continuous_timeseries.budget_compatible_pathways`
+"""
+
+from __future__ import annotations
+
+import sys
+from contextlib import nullcontext as does_not_raise
+from unittest.mock import patch
+
+import numpy as np
+import openscm_units
+import pint.testing
+import pytest
+
+from continuous_timeseries import InterpolationOption, Timeseries
+from continuous_timeseries import budget_compatible_pathways as ct_bcp
+from continuous_timeseries.exceptions import MissingOptionalDependencyError
+
+UR = openscm_units.unit_registry
+Q = UR.Quantity
+
+
+budget_cases = pytest.mark.parametrize(
+    "budget, budget_start_time, emissions_start",
+    (
+        pytest.param(
+            Q(500, "GtC"),
+            Q(2020, "yr"),
+            Q(10.2, "GtC / yr"),
+            id="ar6_global_budget_like",
+        ),
+        pytest.param(
+            Q(1.5, "GtCO2"),
+            Q(2025, "yr"),
+            Q(453.2, "MtCO2 / yr"),
+            id="country_small_budget_like",
+        ),
+        pytest.param(
+            Q(15.0, "GtCO2"),
+            Q(2025, "yr"),
+            Q(745.0, "MtCO2 / yr"),
+            id="country_bigger_budget_like",
+        ),
+    ),
+)
+
+
+@budget_cases
+@pytest.mark.parametrize(
+    "name_res",
+    (
+        pytest.param(None, id="default_name_res"),
+        pytest.param("overwritten", id="name_res_supplied"),
+    ),
+)
+def test_budget_compatible_linear_pathway_derivation(
+    name_res,
+    budget,
+    budget_start_time,
+    emissions_start,
+):
+    kwargs = {}
+    if name_res is not None:
+        kwargs["name_res"] = name_res
+
+    res = ct_bcp.derive_linear_path(
+        budget=budget,
+        budget_start_time=budget_start_time,
+        emissions_start=emissions_start,
+        **kwargs,
+    )
+
+    if name_res is None:
+        assert res.name == (
+            "Linear emissions\n"
+            f"compatible with budget of {budget:.2f}\n"
+            f"from {budget_start_time:.2f}"
+        )
+
+    else:
+        assert res.name == name_res
+
+    # First value should be input emissions
+    discrete_values = res.discrete.values_at_bounds.values
+    pint.testing.assert_equal(discrete_values[0], emissions_start)
+
+    # Last two values should be zero
+    pint.testing.assert_equal(discrete_values[-2:], 0.0 * emissions_start)
+
+    # Derivative should be constant (i.e. function is linear),
+    # except for the last two timesteps,
+    # which should both have a gradient of zero.
+    derivative_values = (
+        ct_bcp.convert_to_annual_time_axis(res)
+        .differentiate()
+        .discrete.values_at_bounds.values
+    )
+    pint.testing.assert_equal(derivative_values[:-2], derivative_values[0])
+    pint.testing.assert_equal(derivative_values[-2:], derivative_values[0] * 0.0)
+
+    # Ensure that we budget is respected
+    pint.testing.assert_allclose(
+        res.integrate(
+            Q(0, emissions_start.u * budget_start_time.u)
+        ).discrete.values_at_bounds.values[-1],
+        budget,
+    )
+    # Ensure that extrapolating post the net-zero results in zeros
+    net_zero_year = res.time_axis.bounds[1]
+    res_extrap = res.interpolate(
+        np.hstack(
+            [net_zero_year, net_zero_year + 3 * (net_zero_year - budget_start_time)]
+        ),
+        allow_extrapolation=True,
+    )
+    pint.testing.assert_equal(
+        res_extrap.discrete.values_at_bounds.values, 0.0 * emissions_start
+    )
+
+
+@budget_cases
+@pytest.mark.parametrize(
+    "name_res",
+    (
+        pytest.param(None, id="default_name_res"),
+        pytest.param("overwritten", id="name_res_supplied"),
+    ),
+)
+def test_budget_compatible_symmetric_quadratic_pathway_derivation(
+    name_res,
+    budget,
+    budget_start_time,
+    emissions_start,
+):
+    kwargs = {}
+    if name_res is not None:
+        kwargs["name_res"] = name_res
+
+    res = ct_bcp.derive_symmetric_quadratic_path(
+        budget=budget,
+        budget_start_time=budget_start_time,
+        emissions_start=emissions_start,
+        **kwargs,
+    )
+
+    if name_res is None:
+        assert res.name == (
+            "Symmetric quadratic emissions\n"
+            f"compatible with budget of {budget:.2f}\n"
+            f"from {budget_start_time:.2f}"
+        )
+
+    else:
+        assert res.name == name_res
+
+    # First value should be input emissions
+    discrete_values = res.discrete.values_at_bounds.values
+    pint.testing.assert_equal(discrete_values[0], emissions_start)
+
+    # Last two values should be zero
+    pint.testing.assert_equal(discrete_values[-2:], 0.0 * emissions_start)
+
+    # Derivative should start at zero,
+    # then climb to its peak,
+    # then be back at zero by the net-zero year (two quadratics).
+    net_zero_year = res.time_axis.bounds[-2]
+    derivative_times = np.hstack(
+        [
+            budget_start_time,
+            (budget_start_time + net_zero_year) / 2.0,
+            net_zero_year,
+        ]
+    )
+    derivative_values = (
+        res.differentiate()
+        .interpolate(derivative_times)
+        .discrete.values_at_bounds.values
+    )
+    # Zero at start and end
+    pint.testing.assert_equal(0.0 * derivative_values[0], derivative_values[[0, -1]])
+    # Peak gradient
+    nzd = (net_zero_year - budget_start_time) / 2.0
+    exp_peak_gradient = -emissions_start / nzd
+    pint.testing.assert_allclose(derivative_values.min(), exp_peak_gradient)
+
+    # Two symmetric quadratics.
+    exp_a = emissions_start / (2.0 * nzd**2)
+    exp_second_derivative_abs = 2.0 * exp_a
+    second_derivative_times = (
+        np.union1d(
+            derivative_times.to("yr").m,
+            (
+                np.linspace(
+                    budget_start_time.to("yr").m,
+                    net_zero_year.to("yr").m,
+                    100,
+                )
+            ),
+        )
+        * budget_start_time.to("yr").u
+    )
+    second_derivative_values = (
+        res.differentiate()
+        .differentiate()
+        .interpolate(second_derivative_times)
+        .discrete.values_at_bounds.values
+    )
+    pre_mid_point = second_derivative_times < (budget_start_time + net_zero_year) / 2.0
+    pint.testing.assert_allclose(
+        second_derivative_values[pre_mid_point], -exp_second_derivative_abs
+    )
+    pint.testing.assert_allclose(
+        second_derivative_values[~pre_mid_point][:-1], exp_second_derivative_abs
+    )
+    # At net-zero year, second derivative is zero
+    pint.testing.assert_allclose(
+        second_derivative_values[-1], 0.0 * exp_second_derivative_abs
+    )
+
+    # Ensure that we budget is respected
+    pint.testing.assert_allclose(
+        res.integrate(
+            Q(0, emissions_start.u * budget_start_time.u)
+        ).discrete.values_at_bounds.values[-1],
+        budget,
+    )
+    # Ensure that extrapolating post the net-zero results in zeros
+    res_extrap = res.interpolate(
+        np.hstack(
+            [net_zero_year, net_zero_year + 3 * (net_zero_year - budget_start_time)]
+        ),
+        allow_extrapolation=True,
+    )
+    pint.testing.assert_equal(
+        res_extrap.discrete.values_at_bounds.values, 0.0 * emissions_start
+    )
+
+
+# ensure convert_to_annual_constant_emissions works
+# for integer and non-integer start year
+# (integer start, just get steps)
+# ()
+@pytest.mark.parametrize(
+    "start_time_axis, exp_out_time_axis",
+    (
+        pytest.param(
+            Q([2010.0, 2029.5, 2030.0], "yr"),
+            Q(np.arange(2010, 2030 + 1), "yr"),
+            id="integer_start_year_integer_end_year",
+        ),
+        pytest.param(
+            Q([2010.5, 2050.3, 2052.0], "yr"),
+            Q(
+                np.hstack([2010.5, np.arange(2011.0, 2052 + 0.1)]),
+                "yr",
+            ),
+            id="non_integer_start_year_integer_end_year",
+        ),
+        pytest.param(
+            Q([2010.0, 2050.3, 2052.5], "yr"),
+            Q(np.arange(2010, 2053 + 1), "yr"),
+            id="integer_start_year_non_integer_end_year",
+        ),
+        pytest.param(
+            Q([2010.5, 2050.0, 2050.5], "yr"),
+            Q(
+                np.hstack([2010.5, np.arange(2011.0, 2051 + 0.1)]),
+                "yr",
+            ),
+            id="non_integer_start_year_non_integer_end_year",
+        ),
+    ),
+)
+@pytest.mark.parametrize(
+    "name_res",
+    (
+        pytest.param(None, id="default_name_res"),
+        pytest.param("overwritten", id="name_res_supplied"),
+    ),
+)
+def test_convert_to_annual_constant_emissions(
+    name_res, start_time_axis, exp_out_time_axis
+):
+    name_start = "name_start"
+
+    kwargs = {}
+    if name_res is not None:
+        kwargs["name_res"] = name_res
+
+    start = Timeseries.from_arrays(
+        time_axis_bounds=start_time_axis,
+        values_at_bounds=Q([1.0, 0.0, 0.0], "MtCO2 / yr"),
+        interpolation=InterpolationOption.Linear,
+        name=name_start,
+    )
+
+    res = ct_bcp.convert_to_annual_constant_emissions(start, **kwargs)
+
+    if name_res is not None:
+        assert res.name == name_res
+    else:
+        assert res.name == f"{start.name}_annualised"
+
+    pint.testing.assert_equal(res.time_axis.bounds, exp_out_time_axis)
+
+    # Result should be piecewise constant, hence
+    pint.testing.assert_allclose(
+        res.differentiate().discrete.values_at_bounds.values, Q(0.0, "MtCO2 / yr / yr")
+    )
+
+    # Check integral preserved
+    exp_integral = start.integrate(Q(0.0, "MtCO2")).discrete.values_at_bounds.values[-1]
+    pint.testing.assert_allclose(
+        exp_integral,
+        res.integrate(Q(0.0, "MtCO2")).discrete.values_at_bounds.values[-1],
+    )
+    if res.time_axis.bounds.m[0] == int(res.time_axis.bounds.m[0]):
+        # Make sure that simply adding up values gives back our budget
+        np.testing.assert_allclose(
+            exp_integral,
+            np.cumsum(res.discrete.values_at_bounds.values)
+            .to(exp_integral.u / res.time_axis.bounds.u)
+            .m[-1],
+            rtol=1e-3,
+        )
+
+    # else:
+    #     # In this case, the awkward half year means that simple addition doesn't work.
+    #     # There is nothing we can do about that.
+
+
+@pytest.mark.parametrize(
+    "sys_modules_patch, expectation",
+    (
+        pytest.param({}, does_not_raise(), id="scipy_available"),
+        pytest.param(
+            {"scipy": None},
+            pytest.raises(
+                MissingOptionalDependencyError,
+                match=(
+                    "`derive_symmetric_quadratic_path` "
+                    "requires scipy to be installed"
+                ),
+            ),
+            id="scipy_not_available",
+        ),
+    ),
+)
+def test_no_scipy_quadratic(sys_modules_patch, expectation):
+    with patch.dict(sys.modules, sys_modules_patch):
+        with expectation:
+            ct_bcp.derive_symmetric_quadratic_path(
+                budget=Q(450.0, "GtC"),
+                budget_start_time=Q(2010.0, "yr"),
+                emissions_start=Q(2.3, "GtC / yr"),
+            )
diff --git a/tests/integration/test_timeseries_integration.py b/tests/integration/test_timeseries_integration.py
index 03df0ea..ef149db 100644
--- a/tests/integration/test_timeseries_integration.py
+++ b/tests/integration/test_timeseries_integration.py
@@ -513,7 +513,13 @@ def test_discrete(operations_test_case):
 
 
 @operations_test_cases
-@pytest.mark.parametrize("name_res", (None, "overwritten"))
+@pytest.mark.parametrize(
+    "name_res",
+    (
+        pytest.param(None, id="default_name_res"),
+        pytest.param("overwritten", id="name_res_supplied"),
+    ),
+)
 def test_differentiate(operations_test_case, name_res):
     kwargs = {}
     if name_res is not None:
@@ -544,7 +550,13 @@ def test_differentiate(operations_test_case, name_res):
 
 
 @operations_test_cases
-@pytest.mark.parametrize("name_res", (None, "overwritten"))
+@pytest.mark.parametrize(
+    "name_res",
+    (
+        pytest.param(None, id="default_name_res"),
+        pytest.param("overwritten", id="name_res_supplied"),
+    ),
+)
 def test_integrate(operations_test_case, name_res):
     kwargs = {}
     if name_res is not None:
@@ -840,7 +852,13 @@ def get_test_update_interpolation_piecewise_constant_to_higher_order_cases() ->
         ),
     ),
 )
-@pytest.mark.parametrize("name_res", (None, "overwritten"))
+@pytest.mark.parametrize(
+    "name_res",
+    (
+        pytest.param(None, id="default_name_res"),
+        pytest.param("overwritten", id="name_res_supplied"),
+    ),
+)
 def test_update_interpolation(  # noqa: PLR0913
     name_res,
     start_interp,
@@ -955,7 +973,13 @@ def get_test_update_interpolation_integral_preserving_cases() -> tuple[pytest.pa
     "start_interp, end_interp, exp_successful, kwargs",
     get_test_update_interpolation_integral_preserving_cases(),
 )
-@pytest.mark.parametrize("name_res", (None, "overwritten"))
+@pytest.mark.parametrize(
+    "name_res",
+    (
+        pytest.param(None, id="default_name_res"),
+        pytest.param("overwritten", id="name_res_supplied"),
+    ),
+)
 def test_update_interpolation_integral_preserving(
     name_res, start_interp, end_interp, exp_successful, kwargs
 ):
diff --git a/uv.lock b/uv.lock
index b558697..cc5011d 100644
--- a/uv.lock
+++ b/uv.lock
@@ -336,7 +336,7 @@ wheels = [
 
 [[package]]
 name = "continuous-timeseries"
-version = "0.2.0"
+version = "0.2.1a1"
 source = { editable = "." }
 dependencies = [
     { name = "attrs" },
@@ -427,6 +427,7 @@ docs = [
 tests = [
     { name = "coverage" },
     { name = "ipython" },
+    { name = "openscm-units" },
     { name = "pytest" },
     { name = "pytest-cov" },
     { name = "pytest-regressions" },
@@ -511,6 +512,7 @@ docs = [
 tests = [
     { name = "coverage", specifier = ">=7.6.0" },
     { name = "ipython", specifier = ">=8.18.1" },
+    { name = "openscm-units", specifier = ">=0.6.3" },
     { name = "pytest", specifier = ">=8.3.3" },
     { name = "pytest-cov", specifier = ">=5.0.0" },
     { name = "pytest-regressions", specifier = ">=2.6.0" },