diff --git a/docs/methods/properties/bandstructure_gw.md b/docs/methods/properties/bandstructure_gw.md
index 4612dc1290..81ff7b4a4b 100644
--- a/docs/methods/properties/bandstructure_gw.md
+++ b/docs/methods/properties/bandstructure_gw.md
@@ -5,8 +5,8 @@ GW for molecules/condensed phase systems with CP2K. In DFT, the energy of a mole
 corresponds to an eigenvalue of the Kohn-Sham matrix. In GW, the procedure for getting the level
 energies is to first perform a DFT calculation (commonly with the PBE or PBE0 functional) to get the
 molecular orbital wavefunctions and then compute a new GW energy for the molecular orbitals of
-interest. For an introduction into the concept of GW, please have a look at [Golze2019](#Golze2019)
-or read Sec. II and the introduction to Sec. III in [Hueser2013](#Hueser2013).
+interest. For an introduction into the concept of GW, please have a look at [](#Golze2019) or read
+Sec. II and the introduction to Sec. III in [](#Hueser2013).
 
 The GW implementation in CP2K is based on the developments described in [](#Wilhelm2016). The
 computational cost of GW is comparable to RPA and MP2 total energy calculations and therefore high.
@@ -27,14 +27,14 @@ calculation is used for computing the molecular orbitals which can be seen from
 [XC_FUNCTIONAL PBE](#CP2K_INPUT.FORCE_EVAL.DFT.XC.XC_FUNCTIONAL.SECTION_PARAMETERS). The input
 parameters for G0W0 are commented below. While the calculation is running, you can look up the
 G0W0@PBE value for the highest occupied molecular orbital (HOMO) and the lowest unoccupied molecular
-orbital (LUMO) in [vanSetten2015](#vanSetten2015) \[Table 2 and 3 (column AIMS-P16/TM-no RI,
-molecule index 76) which should be -11.97 eV and 2.37 eV for HOMO and LUMO, respectively\]. CP2K
-should be able to exactly reproduce these values. In the output of CP2K, the G0W0@PBE results are
-listed after the SCF after the headline *GW quasiparticle energies*.
+orbital (LUMO) in [](#vanSetten2015) \[Table 2 and 3 (column AIMS-P16/TM-no RI, molecule index 76)
+which should be -11.97 eV and 2.37 eV for HOMO and LUMO, respectively\]. CP2K should be able to
+exactly reproduce these values. In the output of CP2K, the G0W0@PBE results are listed after the SCF
+after the headline *GW quasiparticle energies*.
 
 For checking the basis set convergence, we refer to a detailed analysis in [](#Wilhelm2016) in Fig.
 2 for benzene. An extensive table of basis set extrapolated GW levels can be found in
-[vanSetten2015](#vanSetten2015) (column EXTRA).
+[](#vanSetten2015) (column EXTRA).
 
 ```none
 &FORCE_EVAL
@@ -169,19 +169,19 @@ with 1/N_basis   -11.97 +/- 0.02       1.71 +/- 0.29
 GW100            -12.05                2.01
 ```
 
-For the extrapolation, two schemes have been used as described in [vanSetten2015](#vanSetten2015)
-and its supporting information. The first scheme employs a linear fit on the HOMO or LUMO values
-when they are plotted against the inverse cardinal number $N_\text{card}$ of the basis set while the
-second scheme extrapolates versus the inverse number of basis functions $N_\text{basis}$ which can
-be computed as sum of the number of occupied orbitals and the number of virtual orbitals as printed
-in RI_INFO in the output. You can check the extrapolation from the table above with your tool of
+For the extrapolation, two schemes have been used as described in [](#vanSetten2015) and its
+supporting information. The first scheme employs a linear fit on the HOMO or LUMO values when they
+are plotted against the inverse cardinal number $N_\text{card}$ of the basis set while the second
+scheme extrapolates versus the inverse number of basis functions $N_\text{basis}$ which can be
+computed as sum of the number of occupied orbitals and the number of virtual orbitals as printed in
+RI_INFO in the output. You can check the extrapolation from the table above with your tool of
 choice.
 
 The basis set extrapolated values from the table above deviate from the values reported in
 [vanSetten2015](#vanSetten2015), probably because only two basis sets (def2-TZVP, def2-QZVP) have
-been used in [vanSetten2015](#vanSetten2015) for the extrapolation. The extrapolation for the LUMO
-is not working well because one would need much more diffuse functions to represent unbound
-electronic levels (with positive energy).
+been used in [](#vanSetten2015) for the extrapolation. The extrapolation for the LUMO is not working
+well because one would need much more diffuse functions to represent unbound electronic levels (with
+positive energy).
 
 Often, the HOMO-LUMO gap is of interest. In this case, augmented basis sets (e.g. from the EMSL
 database) can offer an alternative for very fast basis set convergence, see also Fig. 2b in
@@ -204,9 +204,17 @@ The G0W0@PBE HOMO value of the H2O molecule (~ -12.0 eV) is not in good agreemen
 experimental ionization potential (12.62 eV). Benchmarks on molecules and solids indicate that
 self-consistency of eigenvalues in the Green's function G improves the agreement between the GW
 calculation and experiment. This scheme is called GW0@PBE and is our favorite GW flavor so far
-(experience based on nano-sized aromatic molecules). You can also check
-[this video](https://www.youtube.com/watch?v=1vUuethWhbs&t=5563s) or [Golze2019](#Golze2019) on
-which DFT starting functional and which self-consistency scheme could be good for your system.
+(experience based on nano-sized aromatic molecules). You can also check the following video or
+[](#Golze2019) on which DFT starting functional and which self-consistency scheme could be good for
+your system.
+
+```{youtube} 1vUuethWhbs
+---
+url_parameters: ?start=5563
+align: center
+privacy_mode:
+---
+```
 
 You can run GW0 calculations in CP2K by putting
 
diff --git a/docs/methods/properties/optical/bethe-salpeter.md b/docs/methods/properties/optical/bethe-salpeter.md
index c75d26bcc8..59424d0db7 100644
--- a/docs/methods/properties/optical/bethe-salpeter.md
+++ b/docs/methods/properties/optical/bethe-salpeter.md
@@ -2,19 +2,18 @@
 
 In this section, we discuss the basics for computing optical properties of molecules using the
 Bethe-Salpeter equation (BSE) in CP2K. The BSE enables the computation of electronic excitation
-energies and optical absorption spectra, for a review, see \[[Blase2018](#Blase2018),
-[Blase2020](#Blase2020), [Bruneval2015](#Bruneval2015), [Sander2015](#Sander2015)\]. In this howto,
-we describe in Sec. [1](#header-theory) the theory and implementation of BSE, in Sec.
-[2](#header-input) the BSE input keywords and in Sec. [3](#header-example) a full CP2K input file of
-a BSE calculation and the corresponding output.
+energies and optical absorption spectra, for a review, see \[[](#Blase2018), [](#Blase2020),
+[](#Bruneval2015), [](#Sander2015)\]. In this howto, we describe in Sec. [1](#header-theory) the
+theory and implementation of BSE, in Sec. [2](#header-input) the BSE input keywords and in Sec.
+[3](#header-example) a full CP2K input file of a BSE calculation and the corresponding output.
 
 (header-theory)=
 
 ## 1. Theory and implementation of BSE
 
 A central goal of a BSE calculation is to compute electronic excitation energies
-$\Omega^{(n)}, n=1,2,\ldots$ (cf. Refs. \[[Blase2018](#Blase2018), [Blase2020](#Blase2020)\] for
-more usecases and details).
+$\Omega^{(n)}, n=1,2,\ldots$ (cf. Refs. \[[](#Blase2018), [](#Blase2020)\] for more usecases and
+details).
 
 The following ingredients are necessary for computing $\Omega^{(n)}$:
 
@@ -24,7 +23,7 @@ The following ingredients are necessary for computing $\Omega^{(n)}$:
 - $GW$ eigenvalues $\varepsilon_i^{GW}$ and $\varepsilon_a^{GW}$ of corresponding KS orbitals.
 
 In CP2K, it is possible to use $G_0W_0$, ev$GW_0$ or ev$GW$ eigenvalues, see details in [GW] and in
-Ref. \[[Golze2019](#Golze2019)\], i.e. we perform BSE@$G_0W_0$/ev$GW_0$/ev$GW$@DFT. Thus, also input
+Ref. \[[](#Golze2019)\], i.e. we perform BSE@$G_0W_0$/ev$GW_0$/ev$GW$@DFT. Thus, also input
 parameters for a DFT and $GW$ calculation can be given (see full input in Sec.
 [3.1](#header-input-file)). We obtain optical properties from BSE solving the following generalized
 eigenvalue problem that involves the block matrix $ABBA$:
@@ -33,7 +32,7 @@ $$\left( \begin{array}{cc}A &  B\\B &  A\end{array} \right)\left( \begin{array}{
 
 We abbreviate $A$ and $B$ as matrices with index $A_{ia,jb}$, i.e. they have
 $N_\mathrm{occ}N_\mathrm{empty}$ rows and $N_\mathrm{occ}N_\mathrm{empty}$ columns. The entries of
-$A$ and $B$ are given by \[[Blase2018](#Blase2018)\]
+$A$ and $B$ are given by \[[](#Blase2018)\]
 
 $$ \begin{align}
     A_{ia,jb} &= (\varepsilon_a^{GW}-\varepsilon_i^{GW})\delta_{ij}\delta_{ab} + \alpha^\mathrm{S/T}
@@ -47,7 +46,7 @@ interaction and $W_{pq,rs}(\omega=0)$ the statically ($\omega=0$) screened Coulo
 where $p,q,r,s \in [ 1, N_\mathrm{occ}+N_\mathrm{empty}]$ are KS orbital indices.
 $(\mathbf{X}^{(n)},\mathbf{Y}^{(n)})$ with elements $X_{ia}^{(n)}$ and $Y_{ia}^{(n)}$ are the
 eigenvectors of the excitation $n$ which relate to the wavefunction of the electronic excitation
-\[[Blase2020](#Blase2020)\],
+\[[](#Blase2020)\],
 
 $$ \begin{align}
 \Psi_\text{excitation}^{(n)}(\mathbf{r}_e,\mathbf{r}_h) = \sum_{ia} X_{ia}^{(n)} \varphi_i(\mathbf{r}_h) \varphi_a(\mathbf{r}_e) + Y_{ia}^{(n)} \varphi_i(\mathbf{r}_e) \varphi_a(\mathbf{r}_h) \quad ,
@@ -114,7 +113,7 @@ In the upper GW/BSE section, the following keywords have been used:
   Setting a small `ENERGY_CUTOFF_OCC` drastically reduces the computation time and the memory
   consumption, but also might affect the computed excitation energies $\Omega^{(n)}$. Recommended to
   use for large systems with more than 30 atoms, but we recommend a careful convergence test by
-  increasing `ENERGY_CUTOFF_OCC` and observing the effect on $\Omega^{(n)}$ \[[Liu2020](#Liu2020)\].
+  increasing `ENERGY_CUTOFF_OCC` and observing the effect on $\Omega^{(n)}$ \[[](#Liu2020)\].
 
 - [ENERGY_CUTOFF_EMPTY](#CP2K_INPUT.FORCE_EVAL.DFT.XC.WF_CORRELATION.RI_RPA.GW.BSE.ENERGY_CUTOFF_EMPTY)
   $E_\text{cut}^\text{empty}$: Analogous to `ENERGY_CUTOFF_OCC`, but for the empty states, i.e. only
@@ -126,9 +125,9 @@ excitation energies:
 
 - [XC_FUNCTIONAL](#CP2K_INPUT.FORCE_EVAL.DFT.XC.XC_FUNCTIONAL): Choose between one of the available
   xc-functionals. The starting point can have a profound influence on the excitation energies
-  \[[Bruneval2015](#Bruneval2015), [Gui2018](#Gui2018), [Jacquemin2017](#Jacquemin2017)\]. We either
-  recommend the PBE functional as DFT starting point when using BSE@ev$GW_0$@PBE or the PBE0
-  functional, when using BSE@$G_0W_0$@PBE0 (see also
+  \[[Bruneval2015](#Bruneval2015), [](#Gui2018), [](#Jacquemin2017)\]. We either recommend the PBE
+  functional as DFT starting point when using BSE@ev$GW_0$@PBE or the PBE0 functional, when using
+  BSE@$G_0W_0$@PBE0 (see also
   [SELF_CONSISTENCY](#CP2K_INPUT.FORCE_EVAL.DFT.XC.WF_CORRELATION.RI_RPA.GW.SELF_CONSISTENCY)).
 
 - [BASIS_SET](#CP2K_INPUT.FORCE_EVAL.SUBSYS.KIND.BASIS_SET): Specify the basis set, which affects
@@ -144,7 +143,7 @@ $N_\mathrm{empty}$ and the estimated memory consumption from the BSE output. The
 is well parallelized, i.e. you can use several nodes that can provide the memory.
 
 We have benchmarked the numerical precision of our BSE implementation and we found excellent
-agreement within only 10 meV compared to the BSE implementation in FHI aims \[[Liu2020](#Liu2020)\].
+agreement within only 10 meV compared to the BSE implementation in FHI aims \[[](#Liu2020)\].
 
 The current BSE implementation in CP2K works for molecules. The inclusion of periodic boundary
 conditions in a $\Gamma$-only approach and with full $k$-point sampling is work in progress.
@@ -243,7 +242,7 @@ The basis sets `aug-cc-pVDZ` and `aug-cc-pVDZ-RIFIT` in `BASIS-aug` can be obtai
 Set Exchange Library:
 <a href="https://www.basissetexchange.org/basis/aug-cc-pvdz/format/cp2k/?version=1&elements=1" target="_blank">aug-cc-pVDZ</a>,
 <a href="https://www.basissetexchange.org/basis/aug-cc-pvdz-rifit/format/cp2k/?version=1&elements=1" target="_blank">aug-cc-pVDZ-RIFIT</a>.
-The geometry for $\mathrm{H}_2$ was taken from [vanSetten2015](#vanSetten2015).
+The geometry for $\mathrm{H}_2$ was taken from [](#vanSetten2015).
 
 ### 3.2 Output
 
diff --git a/docs/methods/sampling/ehrenfest.md b/docs/methods/sampling/ehrenfest.md
index 4750a3ccd6..a86b6cf34e 100644
--- a/docs/methods/sampling/ehrenfest.md
+++ b/docs/methods/sampling/ehrenfest.md
@@ -516,8 +516,9 @@ excited state for the Oxygen 1s is degenerate: there are two available excited s
 the $xy$-plane. That is why a third projection is requested which uses the other excited state
 proposed by the XAD_TDP module. Therefore, the excited state should be understood as the sum over
 the two equivalent excited states:
-\$$
+
+$$
 n_{Exc}(t) = \rho_{\omega}(t) + \rho_{\omega'}(t)
 $$
-where $\\omega'$ stands for the other equivalent excited state, with $\\omega=\\omega'=529\$
-eV.
+
+where $\omega'$ stands for the other equivalent excited state, with $\omega=\omega'=529$ eV.