diff --git a/docs/notebooks/bose_hubbard.ipynb b/docs/notebooks/bose_hubbard.ipynb
new file mode 100644
index 0000000..808ec29
--- /dev/null
+++ b/docs/notebooks/bose_hubbard.ipynb
@@ -0,0 +1,3353 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "id": "6c56b043",
+ "metadata": {},
+ "source": [
+ "# Bose-Hubbard model"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "22cefce0",
+ "metadata": {},
+ "source": [
+ "This notebook builds a **Bose–Hubbard** Hamiltonian over a graph and minimizes its energy with a **photonic variational circuit** in Optyx. We’ll:\n",
+ " \n",
+ "1) prepare a photonic ansatz, \n",
+ "2) define creation/annihilation and number operators on a single photonic mode, \n",
+ "3) assemble the Hamiltonian \\(H\\) from a NetworkX graph using *function syntax*, \n",
+ "4) evaluate $E(\\boldsymbol\\theta)=\\langle\\psi|H|\\psi\\rangle$ and its gradients, \n",
+ "5) run an optimiser with a decaying learning rate, and \n",
+ "6) plot convergence"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "11e739e8",
+ "metadata": {},
+ "source": [
+ "## Ansatz\n",
+ "Define the variational ansatz:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "id": "12d61180-27de-4d3a-bad6-73a0e94c8b16",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ "\n",
+ "\n",
+ "\n"
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "from optyx import photonic, classical\n",
+ "\n",
+ "circuit = photonic.Create(1, 1, 1) >> photonic.ansatz(3, 4) >> photonic.Id(2) @ classical.Select(1)\n",
+ "circuit.draw()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "685b16c4",
+ "metadata": {},
+ "source": [
+ "## The model"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "a49ddabf",
+ "metadata": {},
+ "source": [
+ "Define the creation/annihilation channels for a single photonic mode; we’ll reuse these to place operators on specific lattice sites.\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "5736ad24",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from optyx import Channel\n",
+ "from optyx.core.diagram import mode\n",
+ "from optyx.core.zw import W, Create, Select\n",
+ "\n",
+ "creation_op = Channel(\n",
+ " \"a†\",\n",
+ " Create(1) @ mode >> W(2).dagger()\n",
+ ")\n",
+ "\n",
+ "annihilation_op = Channel(\n",
+ " \"a\",\n",
+ " W(2) >> Select(1) @ mode\n",
+ ")"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "9e720322",
+ "metadata": {},
+ "source": [
+ "We’ll construct the Bose–Hubbard Hamiltonian $H(t,U,\\mu)$ on an arbitrary graph using function syntax so each term edits only its intended wire(s).\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "77676e33",
+ "metadata": {},
+ "source": [
+ "![Bose Hubbard model](./bose_hubbard.png \"Bose Hubbard model\")"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "1fbf6443",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import networkx as nx\n",
+ "from optyx import Diagram, qmode\n",
+ "from optyx.photonic import NumOp, Scalar\n",
+ "\n",
+ "def bose_hubbard_from_graph(\n",
+ " graph: nx.Graph,\n",
+ " t: float,\n",
+ " mu: float,\n",
+ " U: float\n",
+ "):\n",
+ " nodes = sorted(graph.nodes())\n",
+ " idx = {u: i for i, u in enumerate(nodes)}\n",
+ " N = len(nodes)\n",
+ "\n",
+ " H = None\n",
+ "\n",
+ " # hopping: -t (a_i^dagger a_j + a_j^dagger a_i)\n",
+ " for u, v in graph.edges():\n",
+ " i, j = idx[u], idx[v]\n",
+ "\n",
+ " @Diagram.from_callable(dom=qmode**N, cod=qmode**N)\n",
+ " def hop_ij(*in_wires):\n",
+ " # a_i^dagger a_j\n",
+ " out = list(in_wires)\n",
+ " out[i] = creation_op(out[i])\n",
+ " out[j] = annihilation_op(out[j])\n",
+ " Scalar(-t)()\n",
+ " return tuple(out)\n",
+ "\n",
+ " @Diagram.from_callable(dom=qmode**N, cod=qmode**N)\n",
+ " def hop_ji(*in_wires):\n",
+ " # a_j^dagger a_i\n",
+ " out = list(in_wires)\n",
+ " out[j] = creation_op(out[j])\n",
+ " out[i] = annihilation_op(out[i])\n",
+ " Scalar(-t)()\n",
+ " return tuple(out)\n",
+ "\n",
+ " H = (hop_ij + hop_ji) if H is None else (H + hop_ij + hop_ji)\n",
+ "\n",
+ " # onsite interaction: (U/2) a_i^dagger a_i^dagger a_i a_i\n",
+ " for u in nodes:\n",
+ " i = idx[u]\n",
+ "\n",
+ " @Diagram.from_callable(dom=qmode**N, cod=qmode**N)\n",
+ " def quartic_i(*in_wires, i=i):\n",
+ " out = list(in_wires)\n",
+ " w = out[i]\n",
+ " w = creation_op(w)\n",
+ " w = creation_op(w)\n",
+ " w = annihilation_op(w)\n",
+ " w = annihilation_op(w)\n",
+ " out[i] = w\n",
+ " Scalar(U/2)()\n",
+ " return tuple(out)\n",
+ "\n",
+ " H = quartic_i if H is None else (H + quartic_i)\n",
+ "\n",
+ " # -mu n_i\n",
+ " for u in nodes:\n",
+ " i = idx[u]\n",
+ "\n",
+ " @Diagram.from_callable(dom=qmode**N, cod=qmode**N)\n",
+ " def n_i(*in_wires, i=i):\n",
+ " out = list(in_wires)\n",
+ " out[i] = NumOp()(out[i])\n",
+ " Scalar(-mu)()\n",
+ " return tuple(out)\n",
+ "\n",
+ " H = n_i if H is None else (H + n_i)\n",
+ "\n",
+ " return H\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "770d29fb",
+ "metadata": {},
+ "source": [
+ "Assemble hopping, on-site interaction, and chemical-potential terms into one Diagram representing the full Hamiltonian on $qmode^{\\otimes N}$. Start with a 2-site chain and some parameters:\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 5,
+ "id": "900d0c48",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import networkx as nx\n",
+ "\n",
+ "graph = nx.path_graph(2) # 2 sites\n",
+ "\n",
+ "t, U, mu = 0.10, 4.0, 2.0\n",
+ "\n",
+ "hamiltonian = bose_hubbard_from_graph(graph, t, mu, U)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 6,
+ "id": "bb26b7d4",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ "\n",
+ "\n",
+ "\n"
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "hamiltonian.draw()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 7,
+ "id": "8a02308c-aba4-4a39-9325-faee4c1f0c2d",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "expectation = circuit >> hamiltonian >> circuit.dagger()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "c76edfb7",
+ "metadata": {},
+ "source": [
+ "## Optimisation"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 8,
+ "id": "f4e5462d-464d-4b09-8574-6867522b61ea",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from optyx.core.backends import PermanentBackend\n",
+ "\n",
+ "\n",
+ "def to_float(x):\n",
+ " if isinstance(x, complex):\n",
+ " assert x.imag < 1e-8, x\n",
+ " return x.real\n",
+ " return x\n",
+ "\n",
+ "free_syms = list(expectation.free_symbols)\n",
+ "\n",
+ "f_exp = lambda xs: to_float(\n",
+ " expectation.lambdify(*free_syms)(*xs)\n",
+ " .eval(PermanentBackend())\n",
+ " .tensor\n",
+ " .array\n",
+ ")\n",
+ "\n",
+ "def d_f_exp(xs):\n",
+ " return [\n",
+ " expectation.grad(s).lambdify(*free_syms)(*xs)\n",
+ " .eval(PermanentBackend())\n",
+ " .tensor\n",
+ " .array\n",
+ " for s in free_syms\n",
+ " ]"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "5adb53ee",
+ "metadata": {},
+ "source": [
+ "Run a short gradient-descent loop with a decaying learning rate to drive the energy down."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 38,
+ "id": "24355282-4c9c-40b3-b643-461981f9475f",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "from tqdm import tqdm\n",
+ "\n",
+ "xs = []\n",
+ "fxs = []\n",
+ "dfxs = []\n",
+ "\n",
+ "def optimize(x0):\n",
+ " x = x0\n",
+ " lr = 5\n",
+ " steps = 10\n",
+ " for _ in tqdm(range(steps)):\n",
+ " fx = f_exp(x)\n",
+ " dfx = d_f_exp(x)\n",
+ "\n",
+ " xs.append(x[::])\n",
+ " fxs.append(fx)\n",
+ " dfxs.append(dfx)\n",
+ " for i, dfxx in enumerate(dfx):\n",
+ " x[i] = to_float(x[i] - lr * dfxx)\n",
+ "\n",
+ " lr *= 0.2*(i**(1/6)) # make lr smaller with each step\n",
+ "\n",
+ " xs.append(x[::])\n",
+ " fxs.append(f_exp(x))\n",
+ " dfxs.append(d_f_exp(x))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 39,
+ "id": "cde44850",
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stderr",
+ "output_type": "stream",
+ "text": [
+ "100%|██████████| 10/10 [05:31<00:00, 33.16s/it]\n"
+ ]
+ }
+ ],
+ "source": [
+ "optimize([2]*len(free_syms))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 40,
+ "id": "088c8fa8-76b3-45bd-a5dd-b337be0d33b1",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ "\n",
+ "\n",
+ "\n"
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "import matplotlib.pyplot as plt\n",
+ "import seaborn as sns\n",
+ "sns.set_theme()\n",
+ "\n",
+ "fig, axs = plt.subplots(1, 1, figsize=(10, 7))\n",
+ "\n",
+ "axs.plot(range(len(xs)), fxs, c=\"#0072B2\", marker='o')\n",
+ "axs.set_xlabel('Iteration', fontsize=18)\n",
+ "axs.set_ylabel('Expected Energy', fontsize=18)\n",
+ "axs.grid(True)\n",
+ "axs.tick_params(axis='both', which='major', labelsize=16)\n",
+ "axs.tick_params(axis='both', which='minor', labelsize=16)\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "92dfc684",
+ "metadata": {},
+ "outputs": [],
+ "source": []
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": ".venv",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.12.3"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/docs/notebooks/bose_hubbard.png b/docs/notebooks/bose_hubbard.png
new file mode 100644
index 0000000..ea43a2e
Binary files /dev/null and b/docs/notebooks/bose_hubbard.png differ
diff --git a/docs/notebooks/bosonic-vqe-2.ipynb b/docs/notebooks/bosonic-vqe-2.ipynb
index d028b85..dc860e5 100644
--- a/docs/notebooks/bosonic-vqe-2.ipynb
+++ b/docs/notebooks/bosonic-vqe-2.ipynb
@@ -17,7 +17,7 @@
" \n",
" \n",
" \n",
- " 2025-09-05T20:09:47.097250\n",
+ " 2025-10-09T17:20:26.561086\n",
" image/svg+xml\n",
" \n",
" \n",
@@ -46,62 +46,62 @@
"L 223.2 7.2 \n",
"L 7.2 7.2 \n",
"z\n",
- "\" clip-path=\"url(#pd1650474ac)\" style=\"fill: #ffffff; stroke: #ffffff; stroke-linejoin: miter\"/>\n",
+ "\" clip-path=\"url(#p605fbb3927)\" style=\"fill: #ffffff; stroke: #ffffff; stroke-linejoin: miter\"/>\n",
" \n",
" \n",
" \n",
+ "\" clip-path=\"url(#p605fbb3927)\" style=\"fill: none; stroke: #000000; stroke-linejoin: miter\"/>\n",
" \n",
" \n",
" \n",
+ "\" clip-path=\"url(#p605fbb3927)\" style=\"fill: none; stroke: #000000; stroke-linejoin: miter\"/>\n",
" \n",
" \n",
" \n",
+ "\" clip-path=\"url(#p605fbb3927)\" style=\"fill: none; stroke: #000000; stroke-linejoin: miter\"/>\n",
" \n",
" \n",
" \n",
+ "\" clip-path=\"url(#p605fbb3927)\" style=\"fill: none; stroke: #000000; stroke-linejoin: miter\"/>\n",
" \n",
" \n",
" \n",
+ "\" clip-path=\"url(#p605fbb3927)\" style=\"fill: none; stroke: #000000; stroke-linejoin: miter\"/>\n",
" \n",
" \n",
" \n",
+ "\" clip-path=\"url(#p605fbb3927)\" style=\"fill: none; stroke: #000000; stroke-linejoin: miter\"/>\n",
" \n",
" \n",
" \n",
+ "\" clip-path=\"url(#p605fbb3927)\" style=\"fill: none; stroke: #000000; stroke-linejoin: miter\"/>\n",
" \n",
" \n",
" \n",
+ "\" clip-path=\"url(#p605fbb3927)\" style=\"fill: none; stroke: #000000; stroke-linejoin: miter\"/>\n",
" \n",
" \n",
" \n",
+ "\" clip-path=\"url(#p605fbb3927)\" style=\"fill: none; stroke: #000000; stroke-linejoin: miter\"/>\n",
" \n",
" \n",
" \n",
+ "\" clip-path=\"url(#p605fbb3927)\" style=\"fill: none; stroke: #000000; stroke-linejoin: miter\"/>\n",
" \n",
" \n",
" \n",
+ "\" clip-path=\"url(#p605fbb3927)\" style=\"fill: none; stroke: #000000; stroke-linejoin: miter\"/>\n",
" \n",
" \n",
" \n",
+ "\" clip-path=\"url(#p605fbb3927)\" style=\"fill: #ffffff; stroke: #000000; stroke-linejoin: miter\"/>\n",
" \n",
" \n",
" \n",
+ "\" clip-path=\"url(#p605fbb3927)\" style=\"fill: #ffffff; stroke: #000000; stroke-linejoin: miter\"/>\n",
" \n",
" \n",
" \n",
+ "\" clip-path=\"url(#p605fbb3927)\" style=\"fill: #ffffff; stroke: #000000; stroke-linejoin: miter\"/>\n",
" \n",
" \n",
" \n",
+ "\" clip-path=\"url(#p605fbb3927)\" style=\"fill: #ffffff; stroke: #000000; stroke-linejoin: miter\"/>\n",
" \n",
" \n",
" \n",
+ "\" clip-path=\"url(#p605fbb3927)\" style=\"fill: #ffffff; stroke: #000000; stroke-linejoin: miter\"/>\n",
" \n",
" \n",
" \n",
+ "\" clip-path=\"url(#p605fbb3927)\" style=\"fill: #ffffff; stroke: #000000; stroke-linejoin: miter\"/>\n",
" \n",
" \n",
" \n",
@@ -878,7 +878,7 @@
" \n",
" \n",
" \n",
- " \n",
+ " \n",
" \n",
" \n",
" \n",
@@ -911,12 +911,12 @@
"\n",
"\n",
- "