diff --git a/docs/source/tutorial/graphene_example.md b/docs/source/tutorial/graphene_example.md
index d6cff065f0b4931d93183a97a506370cc20ea810..636f9c56a371219748e0d5121c433455ec5450a4 100644
--- a/docs/source/tutorial/graphene_example.md
+++ b/docs/source/tutorial/graphene_example.md
@@ -110,6 +110,7 @@ For example, we can compute the charge density wave (CDW) order parameter which
 To calculate operator expectation values, we first need to construct the density matrix via the {autolink}`~meanfi.mf.density_matrix` function.
 We then feed it into {autolink}`~meanfi.observables.expectation_value` function together with the operator we want to measure.
 In this case, we compute the CDW order parameter by measuring the expectation value of the $\sigma_z$ operator acting on the graphene sublattice degree of freedom.
+
 ```{code-cell} ipython3
 cdw_operator = {(0, 0): np.kron(sz, np.eye(2))}
 
diff --git a/docs/source/tutorial/hubbard_1d.md b/docs/source/tutorial/hubbard_1d.md
index f32933cad3b3225167ef59598c4961e5d4d8f60d..dfc64d26a99a02e53bf675baecfc94fa9141c562 100644
--- a/docs/source/tutorial/hubbard_1d.md
+++ b/docs/source/tutorial/hubbard_1d.md
@@ -10,6 +10,7 @@ kernelspec:
   language: python
   name: python3
 ---
+
 # 1D Hubbard model
 
 ## Background physics
@@ -48,6 +49,7 @@ import numpy as np
 import matplotlib.pyplot as plt
 import meanfi
 ```
+
 Now let us translate the non-interacting Hamiltonian $\hat{H_0}$ defined above into the basic input format for the package: a **tight-binding dictionary**.
 The tight-binding dictionary is a python dictionary where the keys are tuples of integers representing the hopping vectors and the values are the hopping matrices.
 For example, a key `(0,)` represents the onsite term in one dimension and a key `(1,)` represents the hopping a single unit cell to the right.
@@ -59,6 +61,7 @@ Thus our non-interacting Hamiltonian  becomes:
 hopp = np.kron(np.array([[0, 1], [0, 0]]), np.eye(2))
 h_0 = {(0,): hopp + hopp.T.conj(), (1,): hopp, (-1,): hopp.T.conj()}
 ```
+
 Here `hopp` is the hopping matrix which we define as a kronecker product between sublattice and spin degrees of freedom: `np.array([[0, 1], [0, 0]])` corresponds to the hopping between sublattices and `np.eye(2)` leaves the spin degrees of freedom unchanged.
 In the corresponding tight-binding dictionary `h_0`, the key `(0,)` contains hopping within the unit cell and the keys `(1,)` and `(-1,)` correspond to the hopping between the unit cells to the right and left respectively.
 
@@ -101,6 +104,7 @@ Here `s_x` is the Pauli matrix acting on the spin degrees of freedom, which ensu
 
 To combine the non-interacting and interaction Hamiltonians, we use the {autolink}`~meanfi.model.Model` class.
 In addition to the Hamiltonians, we also need to specify the filling of the system --- the number of electrons per unit cell.
+
 ```{code-cell} ipython3
 filling = 2
 full_model = meanfi.Model(h_0, h_int, filling)
@@ -178,6 +182,7 @@ def compute_phase_diagram(
 
     return np.asarray(gaps, dtype=float)
 
+
 Us = np.linspace(0, 4, 40, endpoint=True)
 gaps = compute_phase_diagram(Us=Us, nk=20, nk_dense=100)