diff --git a/notebooks/layer_fidelity.ipynb b/notebooks/layer_fidelity.ipynb
new file mode 100644
index 0000000..dcf5358
--- /dev/null
+++ b/notebooks/layer_fidelity.ipynb
@@ -0,0 +1,662 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "alleged-prairie",
+   "metadata": {},
+   "source": [
+    "# Run Layer Fidelity (EPLG100) experiment using qiskit-experiments"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "b9b9fbb6",
+   "metadata": {},
+   "source": [
+    "`pip install \"qiskit-experiments>=0.7\" qiskit-ibm-runtime` is required before running this notebook."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "f49e876a",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from qiskit_experiments.library.randomized_benchmarking import LayerFidelity"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "78f23522",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "import rustworkx as rx\n",
+    "import matplotlib.pyplot as plt"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "detailed-limitation",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from qiskit_ibm_runtime import QiskitRuntimeService\n",
+    "service = QiskitRuntimeService()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "fa501105",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "#specify backend and two-qubit gate to be layered \n",
+    "backend = service.backend('ibm_kyoto')\n",
+    "twoq_gate = \"ecr\"\n",
+    "print(f\"Device {backend.name} Loaded with {backend.num_qubits} qubits\")\n",
+    "print(f\"Two Qubit Gate: {twoq_gate}\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "997d7bc9",
+   "metadata": {},
+   "source": [
+    "### Find the best qubit chain based on reported 2q-gate errors\n",
+    "If you want to benchmark EPLG (error per layered gate) on 100 qubits, set `num_qubits_in_chain==100`"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "7bc47055",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "num_qubits_in_chain = 100\n",
+    "coupling_map = backend.target.build_coupling_map(twoq_gate)\n",
+    "G = coupling_map.graph"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "8615f6d8",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def to_edges(path):\n",
+    "    edges = []\n",
+    "    prev_node = None\n",
+    "    for node in path:\n",
+    "        if prev_node is not None:\n",
+    "            if G.has_edge(prev_node, node):\n",
+    "                edges.append((prev_node, node))\n",
+    "            else:\n",
+    "                edges.append((node, prev_node))\n",
+    "        prev_node = node\n",
+    "    return edges\n",
+    "    \n",
+    "def path_fidelity(path, correct_by_duration: bool=True, readout_scale: float=None):\n",
+    "    \"\"\"Compute an estimate of the total fidelity of 2-qubit gates on a path.\n",
+    "    If `correct_by_duration` is true, each gate fidelity is worsen by\n",
+    "    scale = max_duration / duration, i.e. gate_fidelity^scale.\n",
+    "    If `readout_scale` > 0 is supplied, readout_fidelity^(1/readout_scale)\n",
+    "    for each qubit on the path is multiplied to the total fielity.\n",
+    "    The path is given in node indices form, e.g. [0, 1, 2].\n",
+    "    An external function `to_edges` is used to obtain edge list, e.g. [(0, 1), (1, 2)].\"\"\"\n",
+    "    path_edges = to_edges(path)\n",
+    "    max_duration = max(backend.target[twoq_gate][qs].duration for qs in path_edges)\n",
+    "    \n",
+    "    def gate_fidelity(qpair):\n",
+    "        duration = backend.target[twoq_gate][qpair].duration\n",
+    "        scale = max_duration / duration if correct_by_duration else 1.0\n",
+    "        # 1.25 = (d+1)/d) with d = 4\n",
+    "        return max(0.25, 1 - (1.25 * backend.target[twoq_gate][qpair].error)) ** scale\n",
+    "    \n",
+    "    def readout_fidelity(qubit):\n",
+    "        return max(0.25, 1 - backend.target[\"measure\"][qubit].error)\n",
+    "  \n",
+    "    total_fidelity = np.prod([gate_fidelity(qs) for qs in path_edges])\n",
+    "    if readout_scale:\n",
+    "        total_fidelity *= np.prod([readout_fidelity(q) for q in path]) ** (1/readout_scale)\n",
+    "    return total_fidelity\n",
+    "\n",
+    "def flatten(paths, cutoff=None):  # cutoff not to make run time too large\n",
+    "    return [val for s, ps in paths.items() for t, vals in ps.items() for val in vals[:cutoff] if s < t]\n",
+    "\n",
+    "paths = rx.all_pairs_all_simple_paths(G.to_undirected(multigraph=False), min_depth=num_qubits_in_chain, cutoff=num_qubits_in_chain)\n",
+    "paths = flatten(paths, cutoff=400)\n",
+    "if not paths:\n",
+    "    raise Exception(f\"No qubit chain with length={num_qubits_in_chain} exists in {backend.name}. Try smaller num_qubits_in_chain.\")\n",
+    "\n",
+    "print(f\"Selecting the best from {len(paths)} candidate paths (will take a few minutes)\")\n",
+    "best_qubit_chain = max(paths, key=path_fidelity)\n",
+    "assert(len(best_qubit_chain) == num_qubits_in_chain)\n",
+    "print(f\"Predicted LF from reported 2q-gate EPGs: {path_fidelity(best_qubit_chain)}\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "efc3c14d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "np.array(best_qubit_chain)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "e89c690a",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# decompose the chain into two disjoint layers\n",
+    "all_pairs = to_edges(best_qubit_chain)\n",
+    "two_disjoint_layers = [all_pairs[0::2], all_pairs[1::2]]\n",
+    "two_disjoint_layers"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "inclusive-vietnam",
+   "metadata": {},
+   "source": [
+    "## Define a LayerFidelity experiment"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "rational-branch",
+   "metadata": {},
+   "source": [
+    "- `physical_qubits`: Physical qubits to use in this benchmarking\n",
+    "- `two_qubit_layers`: List of list of qubit pairs (a list of qubit pairs corresponds with a two-qubit gate layer), qubit direction matters\n",
+    "- `lengths`: List of layer lengths (the number of depth points)\n",
+    "- `backend`: Backend to be benchmarked\n",
+    "- `num_samples`: How many random circuits should be sampled for each layer length\n",
+    "- `seed`: Seed for the random number generator\n",
+    "\n",
+    "[Optional]\n",
+    "- `two_qubit_gate`: Two-qubit gate name to use in layers (If not specified, one of 2q-gates supported in the backend is automatically set.)\n",
+    "- `one_qubit_basis_gates`: One-qubit gate names to use for random 1q-Clifford layers (If not specified, all 1q-gates supported in the backend are automatically set.)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "russian-billion",
+   "metadata": {
+    "scrolled": false
+   },
+   "outputs": [],
+   "source": [
+    "%%time\n",
+    "lfexp = LayerFidelity(\n",
+    "    physical_qubits=best_qubit_chain,\n",
+    "    two_qubit_layers=two_disjoint_layers,\n",
+    "    lengths=[2, 4, 8, 16, 30, 50, 70, 100, 150, 200, 300, 500],\n",
+    "    backend=backend,\n",
+    "    num_samples=12,\n",
+    "    seed=42,\n",
+    "    # two_qubit_gate=\"ecr\",\n",
+    "    # one_qubit_basis_gates=[\"rz\", \"sx\", \"x\"],\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "bc9afcb9",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# set maximum number of circuits per job to avoid errors due to too large payload\n",
+    "lfexp.experiment_options.max_circuits = 144"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "50e33867",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "lfexp.experiment_options"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "c033c62c",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "print(f\"Two Qubit Gate: {lfexp.experiment_options.two_qubit_gate}\")\n",
+    "print(f\"One Qubit Gate Set: {lfexp.experiment_options.one_qubit_basis_gates}\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "d19b5a19",
+   "metadata": {},
+   "source": [
+    "### Quick look at a RB circuit\n",
+    "`circuits_geenrator` is useful for that."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "114f5349",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%%time\n",
+    "# look at one of the first three 2Q direct RB circuits quickly\n",
+    "circ_iter = lfexp.circuits_generator()\n",
+    "first_three_circuits = list(next(circ_iter) for _ in range(3))\n",
+    "first_three_circuits[1].draw(output=\"mpl\", style=\"clifford\", idle_wires=False, fold=-1)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "5866f7eb",
+   "metadata": {},
+   "source": [
+    "### Generating all RB circuits takes a bit more time"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a4dfc21d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%%time\n",
+    "# generate all circuits to run\n",
+    "circuits = lfexp.circuits()\n",
+    "print(f\"{len(circuits)} circuits are generated.\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "therapeutic-merchant",
+   "metadata": {},
+   "source": [
+    "## Run the Layer Fidelity experiment\n",
+    "\n",
+    "Run = Generate circuits and submit the job to the backend.\n",
+    "\n",
+    "Set `exp_data.auto_save = True` to store the experiment data using qiskit-ibm-experiment service."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "iraqi-dominican",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%%time\n",
+    "#number of shots per job\n",
+    "nshots = 400\n",
+    "\n",
+    "# Run the LF experiment (generate circuits and submit the job)\n",
+    "exp_data = lfexp.run(shots=nshots)\n",
+    "# exp_data.auto_save = True\n",
+    "print(f\"Run experiment: ID={exp_data.experiment_id} with jobs {exp_data.job_ids}]\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "f0104130",
+   "metadata": {},
+   "source": [
+    "## Wait for the result\n",
+    "Retrieve the result after the job is completed and run analysis (fit RB curves and estimate process fidelity for each 1Q/2Q direct RB)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "banner-trader",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%%time\n",
+    "exp_data.block_for_results()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "white-variation",
+   "metadata": {},
+   "source": [
+    "### Print out the result data as a table"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "fabulous-platinum",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Get the result data as a pandas DataFrame\n",
+    "df = exp_data.analysis_results(dataframe=True)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a3d31324",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# LF and LF of each disjoing set\n",
+    "df[(df.name == \"LF\") | (df.name == \"EPLG\") | (df.name == \"SingleLF\")]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "d80948c0",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# 1Q/2Q process fidelities\n",
+    "df[(df.name == \"ProcessFidelity\")].head()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "1ebf0edf",
+   "metadata": {},
+   "source": [
+    "### Display plots (1Q/2Q simultaneous direct RB curves)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "unlike-python",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Plot 1Q/2Q simultaneous direct RB curves\n",
+    "for i in range(0, 5):\n",
+    "    display(exp_data.figure(i))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "151428c3",
+   "metadata": {},
+   "source": [
+    "### Filter qubit pairs with bad process fidelity "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "3207af8e",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "pfdf = df[(df.name == \"ProcessFidelity\")]\n",
+    "pfdf[pfdf.value < 0.8]"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "28f3ee4c",
+   "metadata": {},
+   "source": [
+    "### Invesitigate analysis results labeled \"bad\" quality"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "a3924a1b",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# find bad quality analysis results\n",
+    "pfdf[pfdf.quality==\"bad\"]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "49f56a2d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Display \"bad\" 1Q/2Q direct RB plots\n",
+    "# for qubits in pfdf[pfdf.quality==\"bad\"].qubits:\n",
+    "#     figname = \"DirectRB_Q\" + \"_Q\".join(map(str, qubits))\n",
+    "#     display(exp_data.figure(figname))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "43060a5a",
+   "metadata": {},
+   "source": [
+    "## Plot Layer Fidelity (EPLG) by chain length"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "530a366a",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# fill Process Fidelity values with zeros\n",
+    "pfdf = pfdf.fillna({\"value\": 0})"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "202188f7",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Compute LF by chain length assuming the first layer is full with 2q-gates\n",
+    "lf_sets, lf_qubits = two_disjoint_layers, best_qubit_chain\n",
+    "full_layer = [None] * (len(lf_sets[0]) + len(lf_sets[1]))\n",
+    "full_layer[::2] = lf_sets[0]\n",
+    "full_layer[1::2] = lf_sets[1]\n",
+    "full_layer = [(lf_qubits[0],)] + full_layer + [(lf_qubits[-1],)]\n",
+    "len(full_layer)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "371d7436",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "assert(len(full_layer)== len(lf_qubits) + 1)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "80d4c824",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "pfs = [pfdf.loc[pfdf[pfdf.qubits==qubits].index[0], 'value'] for qubits in full_layer]\n",
+    "pfs = list(map(lambda x: x.n if x != 0 else 0, pfs))\n",
+    "pfs[0] = pfs[0]**2\n",
+    "pfs[-1] = pfs[-1]**2\n",
+    "np.array(pfs)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "330f6c75",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "len(pfs)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "3e1b2500",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "job = service.job(exp_data.job_ids[0])\n",
+    "JOB_DATE = job.creation_date"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "3eb0afe8",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Approximate 1Q RB fidelities at both ends by the square root of 2Q RB fidelity at both ends.\n",
+    "# For example, if we have [(0, 1), (1, 2), (2, 3), (3, 4)] 2Q RB fidelities and if we want to compute a layer fidelity for [1, 2, 3],\n",
+    "# we approximate the 1Q filedities for (1,) and (3,) by the square root of 2Q fidelities of (0, 1) and (3, 4).\n",
+    "chain_lens = list(range(4, len(pfs), 2))\n",
+    "chain_fids = []\n",
+    "for length in chain_lens:\n",
+    "    w = length + 1  # window size\n",
+    "    fid_w = max(np.sqrt(pfs[s]) * np.prod(pfs[s+1:s+w-1]) * np.sqrt(pfs[s+w-1]) for s in range(len(pfs)-w+1))\n",
+    "    chain_fids.append(fid_w)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "c02610b6",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Plot LF by chain length\n",
+    "plt.title(f\"Backend: {backend.name}, {JOB_DATE.strftime('%Y/%m/%d %H:%M')}\")\n",
+    "plt.plot(\n",
+    "    chain_lens,\n",
+    "    chain_fids,\n",
+    "    marker=\"o\",\n",
+    "    linestyle=\"-\",\n",
+    ")\n",
+    "plt.xlim(0, chain_lens[-1]*1.05)\n",
+    "plt.ylim(0.95*min(chain_fids), 1)\n",
+    "plt.ylabel(\"Layer Fidelity\")\n",
+    "plt.xlabel(\"Chain Length\")\n",
+    "plt.grid()\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "74603316",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Plot EPLG by chain length\n",
+    "num_2q_gates = [length -1 for length in chain_lens]\n",
+    "chain_eplgs = [1 - (fid ** (1 / num_2q)) for num_2q, fid in zip(num_2q_gates, chain_fids)]\n",
+    "plt.title(f\"Backend: {backend.name}, {JOB_DATE.strftime('%Y/%m/%d %H:%M')}\")\n",
+    "plt.plot(\n",
+    "    chain_lens,\n",
+    "    chain_eplgs,\n",
+    "    marker=\"o\",\n",
+    "    linestyle=\"-\",\n",
+    ")\n",
+    "plt.xlim(0, chain_lens[-1]*1.05)\n",
+    "plt.ylabel(\"Error per Layered Gates\")\n",
+    "plt.xlabel(\"Chain Length\")\n",
+    "plt.grid()\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "id": "885f56a7",
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/html": [
+       "<div style='width: 100%; background-color:#d5d9e0;padding-left: 10px; padding-bottom: 10px; padding-right: 10px; padding-top: 5px'><p>&copy; Copyright IBM 2017, 2024.</p><p>This code is licensed under the Apache License, Version 2.0. You may<br>obtain a copy of this license in the LICENSE.txt file in the root directory<br> of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.<p>Any modifications or derivative works of this code must retain this<br>copyright notice, and modified files need to carry a notice indicating<br>that they have been altered from the originals.</p></div>"
+      ],
+      "text/plain": [
+       "<IPython.core.display.HTML object>"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "import datetime\n",
+    "from IPython.display import HTML, display\n",
+    "\n",
+    "def qiskit_copyright(line=\"\", cell=None):\n",
+    "    \"\"\"IBM copyright\"\"\"\n",
+    "    now = datetime.datetime.now()\n",
+    "\n",
+    "    html = \"<div style='width: 100%; background-color:#d5d9e0;\"\n",
+    "    html += \"padding-left: 10px; padding-bottom: 10px; padding-right: 10px; padding-top: 5px'>\"\n",
+    "    html += \"<p>&copy; Copyright IBM 2017, %s.</p>\" % now.year\n",
+    "    html += \"<p>This code is licensed under the Apache License, Version 2.0. You may<br>\"\n",
+    "    html += \"obtain a copy of this license in the LICENSE.txt file in the root directory<br> \"\n",
+    "    html += \"of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.\"\n",
+    "\n",
+    "    html += \"<p>Any modifications or derivative works of this code must retain this<br>\"\n",
+    "    html += \"copyright notice, and modified files need to carry a notice indicating<br>\"\n",
+    "    html += \"that they have been altered from the originals.</p>\"\n",
+    "    html += \"</div>\"\n",
+    "    return display(HTML(html))\n",
+    "\n",
+    "qiskit_copyright()"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "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.10.13"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}