diff --git a/.github/workflows/website.yaml b/.github/workflows/website.yaml
index 7e6c8b71..fb6be00f 100644
--- a/.github/workflows/website.yaml
+++ b/.github/workflows/website.yaml
@@ -22,8 +22,8 @@ jobs:
       - name: Install pip
         run: python -m pip install --upgrade pip
 
-      - name: Install Sphinx Read The Docs
-        run: pip install sphinx
+      - name: Install Sphinx
+        run: pip install "sphinx<8"
 
       - name: Install Read the Docs theme for Sphinx
         run: pip install sphinx-rtd-theme
diff --git a/doc/introduction.rst b/doc/introduction.rst
index f194af75..42f82335 100644
--- a/doc/introduction.rst
+++ b/doc/introduction.rst
@@ -269,9 +269,7 @@ Some options may take several argument values, which can be combined using *sepa
 Initialization File
 ~~~~~~~~~~~~~~~~~~~
 
-.. index::
-
-   --rcfile, initialization file
+.. index:: --rcfile, initialization file
 
 When Neper is run, it starts by reading commands from an initialization file, :file:`$HOME/.neperrc`, if that file exists.  Another initialization file can be specified with option :data:`--rcfile`, before calling a module:
 
@@ -317,9 +315,7 @@ The environment variable is unset by default.
 Conventions Used in This Manual
 -------------------------------
 
-.. index::
-
-   single: dollar sign $, shell prompt
+.. index:: single: dollar sign $, shell prompt
 
 - A command entered at the terminal is shown like this:
 
diff --git a/doc/neper_m.rst b/doc/neper_m.rst
index 3a8d7e31..e33ff946 100644
--- a/doc/neper_m.rst
+++ b/doc/neper_m.rst
@@ -4,8 +4,7 @@
 .. |rasterscale| replace:: :data:`rasterscale(<fact>)` or :data:`rasterscale(<fact_x>,<fact_y>,<fact_z>)`: scale the number of voxels of a raster tessellation, where :data:`<fact>` is the scaling factor that applies in the three directions, and :data:`<fact_x>`, :data:`<fact_y>` and :data:`<fact_z>` are the scaling factors along :data:`x`, :data:`y` and :data:`z`, respectively
 .. |crop| replace:: :data:`crop(<xmin>,<xmax>,<ymin>,<ymax>,<zmin>,<zmax>)`: crop a region of a raster tessellation, where :data:`<xmin>`, :data:`<xmax>`, :data:`<ymin>`, :data:`<ymax>`, :data:`<zmin>` and :data:`<zmax>` are the minimum and maximum positions along :data:`x`, :data:`y` and :data:`z`, respectively. For 2D raster tessellations, the :data:`z` values can be omitted
 
-.. index::
-   single: -M
+.. index:: single: -M
 
 .. _neper_m:
 
@@ -151,11 +150,9 @@ Meshing Options
 
   **Default value**: :data:`tet` in 3D and :data:`tri` in 2D.
 
-.. index::
-   single: -cl
+.. index:: single: -cl
 
-.. index::
-   single: -rcl
+.. index:: single: -rcl
 
 .. option:: -cl or -rcl <length> (secondary option)
 
@@ -382,9 +379,8 @@ Transformation Options
 
   **Default value**: -.
 
-.. index::
-   single: -sortnode
-   single: -sortelt
+.. index:: single: -sortnode
+.. index:: single: -sortelt
 
 .. option:: -sort{node,elt} <sort_expression> (secondary option)
 
@@ -563,9 +559,8 @@ These advanced options set running conditions for the meshing libraries (2D and
 
   **Default value**: :data:`0.02`.
 
-.. index::
-   single: -mesh2dmaxtime
-   single: -mesh3dmaxtime
+.. index:: single: -mesh2dmaxtime
+.. index:: single: -mesh3dmaxtime
 
 .. option:: -mesh{2,3}dmaxtime <time> (secondary option)
 
@@ -574,9 +569,8 @@ These advanced options set running conditions for the meshing libraries (2D and
 
   **Default value**: :data:`1000`.
 
-.. index::
-   single: -mesh2drmaxtime
-   single: -mesh3drmaxtime
+.. index:: single: -mesh2drmaxtime
+.. index:: single: -mesh3drmaxtime
 
 .. option:: -mesh{2,3}drmaxtime <relative_time> (secondary option)
 
diff --git a/doc/neper_s.rst b/doc/neper_s.rst
index 05d64887..55b3ec09 100644
--- a/doc/neper_s.rst
+++ b/doc/neper_s.rst
@@ -1,5 +1,4 @@
-.. index::
-   single: -S
+.. index:: single: -S
 
 .. _neper_s:
 
diff --git a/doc/neper_t.rst b/doc/neper_t.rst
index c9df1083..69738978 100644
--- a/doc/neper_t.rst
+++ b/doc/neper_t.rst
@@ -4,8 +4,7 @@
 .. |crop| replace:: :data:`crop(<xmin>,<xmax>,<ymin>,<ymax>,<zmin>,<zmax>)`: crop a region of a raster tessellation, where :data:`<xmin>`, :data:`<xmax>`, :data:`<ymin>`, :data:`<ymax>`, :data:`<zmin>` and :data:`<zmax>` are the minimum and maximum positions along :data:`x`, :data:`y` and :data:`z`, respectively. For 2D raster tessellations, the :data:`z` values can be omitted
 .. |rasterscale| replace:: :data:`rasterscale(<fact>)` or :data:`rasterscale(<fact_x>,<fact_y>,<fact_z>)`: : scale the number of voxels of a raster tessellation, where :data:`factor` is the scaling factor that applies in the three directions, and :data:`<fact_x>`, :data:`<fact_y>` and :data:`<fact_z>` are the scaling factors along :data:`x`, :data:`y` and :data:`z`, respectively
 
-.. index::
-   single: -T
+.. index:: single: -T
 
 .. _neper_t:
 
@@ -242,7 +241,7 @@ These options can be used to set the cell morphology.
 
   - :data:`diameq`: the equivalent diameter [#size_diameq]_;
 
-  - :data:`sphericity`: the sphericity, and :data:`1-sphericity`: 1 :math:`-` the sphericity (or :data:`circularity` and :data:`1-circularity`). [#sphericity-circularity]_ [#1-sphericity]_
+  - :data:`sphericity`: the sphericity, and :data:`1-sphericity`: 1 :math:`-` the sphericity (or :data:`circularity` and :data:`1-circularity`). [#sphericity-circularity]_ [#one-sphericity]_
 
     All of :data:`size`, :data:`diameq` and :data:`sphericity` (and their variants) can be defined by :ref:`statistical distributions <statistical_distributions>` or cell by cell.  If the number of cells is defined using option :option:`-n`, the :data:`size` or :data:`diameq` distribution is scaled to get the specified number of cells.  At the opposite, if :data:`-n from_morpho` is used, the number of cells is determined from the :data:`size` or :data:`diameq` distribution.  An interval of possible values can also be provided using :data:`interval(<min>,<max>)`.  Cell-by-cell values can be provided using :data:`file(<file_name>)`, where :data:`<file_name>` contains the cell values.  A unique (numeral) value to be assigned to all cells can also be provided directly.
 
@@ -1214,7 +1213,7 @@ References
 
 .. [#sphericity-circularity] Terms :data:`sphericity` and :data:`circularity` apply to 3D and 2D, respectively, but can be used interchangeably. The sphericity of a polyhedron corresponds to the ratio of the surface area of the sphere of equivalent volume to the surface area of the polyhedron. Similarly, the circularity of a polygon corresponds to the ratio of the perimeter of the circle of equivalent surface area to the perimeter of the polygon.
 
-.. [#1-sphericity] The reason behind the :data:`1-sphericity` (or :data:`1-circularity`) variable is that, for a grain growth microstructure, 1 :math:`-` the sphericity follows a lognormal distribution [#CMAME2018]_.
+.. [#one-sphericity] The reason behind the :data:`1-sphericity` (or :data:`1-circularity`) variable is that, for a grain growth microstructure, 1 :math:`-` the sphericity follows a lognormal distribution [CMAME2018]_.
 
 .. [#endianness] Endianness is both written in the tesr file and tested on the system when reading the tesr file, so that the user normally does not have to care about it (even when transferring files across systems).
 
diff --git a/doc/neper_v.rst b/doc/neper_v.rst
index f12d8c61..9587c76f 100644
--- a/doc/neper_v.rst
+++ b/doc/neper_v.rst
@@ -3,8 +3,7 @@
 .. |scale| replace:: :data:`scale(<fact_x>,<fact_y>,<fact_z>)`: scale by :data:`<fact_x>`, :data:`<fact_y>` and :data:`<fact_z>` along :data:`x`, :data:`y` and :data:`z`, respectively
 .. |data_description| replace:: This can be done either directly, by providing one or several values (e.g. the RGB channel values for color) or indirectly, e.g. using scalar values that are converted in color using a given *color scheme*.  In this case (and for real-type data), a scale image is also generated, and its properties (minimum, maximum and tick values) can be specified.
 
-.. index::
-   single: -V
+.. index:: single: -V
 
 .. _neper_v:
 
@@ -159,8 +158,7 @@ Tessellation Data Loading and Rendering Options
 
 The following option enables the definition of the cell data itself (pole figure space only):
 
-.. index::
-   single: -datacell
+.. index:: single: -datacell
 
 .. option:: -datacell <type>:<value>
 
@@ -187,20 +185,19 @@ The following options enable the definition of the properties (color and size) o
 
 For each entity, all attributes can be set, although the may not apply in certain spaces (see option :option:`-space`). Specifically, :data:`-data{cell,poly,face}rad` do not apply in real space.
 
-.. index::
-   single: -datacellcol
-   single: -datapolycol
-   single: -datafacecol
-   single: -dataedgecol
-   single: -datavercol
-   single: -dataseedcol
-   single: -datavoxcol
-   single: -datavoidvoxcol
-   single: -datacrystaledgecol
-   single: -datacelledgecol
-   single: -datavoxedgecol
-   single: -datacrystalcol
-   single: -data...col
+.. index:: single: -datacellcol
+.. index:: single: -datapolycol
+.. index:: single: -datafacecol
+.. index:: single: -dataedgecol
+.. index:: single: -datavercol
+.. index:: single: -dataseedcol
+.. index:: single: -datavoxcol
+.. index:: single: -datavoidvoxcol
+.. index:: single: -datacrystaledgecol
+.. index:: single: -datacelledgecol
+.. index:: single: -datavoxedgecol
+.. index:: single: -datacrystalcol
+.. index:: single: -data...col
 
 .. option:: -data{cell,poly,face,edge,ver,seed,crystal,vox,voidvox,{crystal,cell,vox}edge}col <color>
 
@@ -244,17 +241,16 @@ For each entity, all attributes can be set, although the may not apply in certai
   **Default value**: :data:`id` for cells and crystals, :data:`white` for faces, :data:`black` for edges and vertices, :data:`gray` for seeds, voxels and void voxels, :data:`black` for voxel edges.
 
 
-.. index::
-   single: -datacellcolscheme
-   single: -datapolycolscheme
-   single: -datafacecolscheme
-   single: -dataedgecolscheme
-   single: -datavercolscheme
-   single: -dataseedcolscheme
-   single: -datavoxcolscheme
-   single: -datavoidvoxcolscheme
-   single: -datavoxedgecolscheme
-   single: -data...colscheme
+.. index:: single: -datacellcolscheme
+.. index:: single: -datapolycolscheme
+.. index:: single: -datafacecolscheme
+.. index:: single: -dataedgecolscheme
+.. index:: single: -datavercolscheme
+.. index:: single: -dataseedcolscheme
+.. index:: single: -datavoxcolscheme
+.. index:: single: -datavoidvoxcolscheme
+.. index:: single: -datavoxedgecolscheme
+.. index:: single: -data...colscheme
 
 .. option:: -data{cell,poly,face,edge,ver,seed,crystal,vox}colscheme <col_scheme>
 
@@ -277,14 +273,13 @@ For each entity, all attributes can be set, although the may not apply in certai
 
   **Default value**: :data:`palette` for integer values, :data:`viridis` for real values, and :data:`rodrigues` for crystal orientations.
 
-.. index::
-   single: -datacelltrs
-   single: -datapolytrs
-   single: -datafacetrs
-   single: -dataedgetrs
-   single: -datavertrs
-   single: -dataseedtrs
-   single: -data...trs
+.. index:: single: -datacelltrs
+.. index:: single: -datapolytrs
+.. index:: single: -datafacetrs
+.. index:: single: -dataedgetrs
+.. index:: single: -datavertrs
+.. index:: single: -dataseedtrs
+.. index:: single: -data...trs
 
 .. option:: -data{cell,poly,face,edge,ver,seed,crystal}trs <transparency>
 
@@ -425,9 +420,8 @@ The following options enable the definition of the properties (color, size, etc.
 
   **Default value**: :data:`palette` for integer values, :data:`viridis` for real values, and :data:`rodrigues` for crystal orientations.
 
-.. index::
-   single: -dataelt2dtrs
-   single: -dataelset2dtrs
+.. index:: single: -dataelt2dtrs
+.. index:: single: -dataelset2dtrs
 
 .. option:: -data{elt2d,elset2d}trs <transparency>
 
@@ -669,13 +663,12 @@ Show Options
 
 The following options apply to the full tessellations, mesh, etc.
 
-.. index::
-   single: -showtess
-   single: -showtesr
-   single: -showmesh
-   single: -showmeshslice
-   single: -showcsys
-   single: -showscale
+.. index:: single: -showtess
+.. index:: single: -showtesr
+.. index:: single: -showmesh
+.. index:: single: -showmeshslice
+.. index:: single: -showcsys
+.. index:: single: -showscale
 
 .. option:: -show{tess,tesr,mesh,meshslice,meshsliceelt1d,csys,scale} <logical>
 
@@ -687,16 +680,15 @@ The following options apply to the full tessellations, mesh, etc.
 
 The following option applies to a tessellation or a raster tessellation.
 
-.. index::
-   single: -showcell
-   single: -showpoly
-   single: -showface
-   single: -showedge
-   single: -showver
-   single: -showseed
-   single: -showcrystal
-   single: -showfaceinter
-   single: -showvox
+.. index:: single: -showcell
+.. index:: single: -showpoly
+.. index:: single: -showface
+.. index:: single: -showedge
+.. index:: single: -showver
+.. index:: single: -showseed
+.. index:: single: -showcrystal
+.. index:: single: -showfaceinter
+.. index:: single: -showvox
 
 .. option:: -show{cell,poly,face,edge,ver,seed,crystal,faceinter,vox} <logical>
 
@@ -713,18 +705,17 @@ The following option applies to a tessellation or a raster tessellation.
 
 The following options apply to the entities of the mesh.  The options apply to 3D elements (:data:`elt3d`) and element sets (:data:`elset3d`), 2D elements (:data:`elt2d`) and element sets (:data:`elset2d`), 1D elements (:data:`elt1d`) and element sets (:data:`elset1d`), 0D elements (:data:`elt0d`) and element sets (:data:`elset0d`), and nodes (:data:`nodes`). Also note that the :data:`elt` and :data:`elset` can be used to denote the highest-dimension elements and elsets.
 
-.. index::
-   single: -showelt
-   single: -showelset
-   single: -showelt0d
-   single: -showelt1d
-   single: -showelt2d
-   single: -showelt3d
-   single: -showelset0d
-   single: -showelset1d
-   single: -showelset2d
-   single: -showelset3d
-   single: -shownode
+.. index:: single: -showelt
+.. index:: single: -showelset
+.. index:: single: -showelt0d
+.. index:: single: -showelt1d
+.. index:: single: -showelt2d
+.. index:: single: -showelt3d
+.. index:: single: -showelset0d
+.. index:: single: -showelset1d
+.. index:: single: -showelset2d
+.. index:: single: -showelset3d
+.. index:: single: -shownode
 
 .. option:: -show{elt,elset,elt{0-3}d,elset{0-3}d,node} <logical>