diff --git a/properties/band_gap.py b/properties/band_gap.py
new file mode 100644
index 00000000..c6922a4f
--- /dev/null
+++ b/properties/band_gap.py
@@ -0,0 +1,159 @@
+#
+# Copyright The NOMAD Authors.
+#
+# This file is part of NOMAD. See https://nomad-lab.eu for further info.
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+#
+
+import numpy as np
+from structlog.stdlib import BoundLogger
+import pint
+from typing import Optional
+
+from nomad.metainfo import Quantity, MEnum, Section, Context
+
+from nomad_simulations.physical_property import PhysicalProperty
+
+
+class ElectronicBandGap(PhysicalProperty):
+    """
+    Energy difference between the highest occupied electronic state and the lowest unoccupied electronic state.
+    """
+
+    # ! implement `iri` and `rank` as part of `m_def = Section()`
+
+    iri = 'http://fairmat-nfdi.eu/taxonomy/ElectronicBandGap'
+
+    type = Quantity(
+        type=MEnum('direct', 'indirect'),
+        description="""
+        Type categorization of the electronic band gap. This quantity is directly related with `momentum_transfer` as by
+        definition, the electronic band gap is `'direct'` for zero momentum transfer (or if `momentum_transfer` is `None`) and `'indirect'`
+        for finite momentum transfer.
+
+        Note: in the case of finite `variables`, this quantity refers to all of the `value` in the array.
+        """,
+    )
+
+    momentum_transfer = Quantity(
+        type=np.float64,
+        shape=[2, 3],
+        description="""
+        If the electronic band gap is `'indirect'`, the reciprocal momentum transfer for which the band gap is defined
+        in units of the `reciprocal_lattice_vectors`. The initial and final momentum 3D vectors are given in the first
+        and second element. Example, the momentum transfer in bulk Si2 happens between the Γ and the (approximately)
+        X points in the Brillouin zone; thus:
+            `momentum_transfer = [[0, 0, 0], [0.5, 0.5, 0]]`.
+
+        Note: this quantity only refers to scalar `value`, not to arrays of `value`.
+        """,
+    )
+
+    spin_channel = Quantity(
+        type=np.int32,
+        description="""
+        Spin channel of the corresponding electronic band gap. It can take values of 0 or 1.
+        """,
+    )
+
+    value = Quantity(
+        type=np.float64,
+        unit='joule',
+        description="""
+        The value of the electronic band gap. This value has to be positive, otherwise it will
+        prop an error and be set to None by the `normalize()` function.
+        """,
+    )
+
+    def __init__(
+        self, m_def: Section = None, m_context: Context = None, **kwargs
+    ) -> None:
+        super().__init__(m_def, m_context, **kwargs)
+        self.name = self.m_def.name
+        self.rank = []
+
+    def validate_values(self, logger: BoundLogger) -> Optional[pint.Quantity]:
+        """
+        Validate the electronic band gap `value` by checking if they are negative and sets them to None if they are.
+
+        Args:
+            logger (BoundLogger): The logger to log messages.
+        """
+        value = self.value.magnitude
+        if not isinstance(self.value.magnitude, np.ndarray):  # for scalars
+            value = np.array(
+                [value]
+            )  # ! check this when talking with Lauri and Theodore
+
+        # Set the value to 0 when it is negative
+        if (value < 0).any():
+            logger.error('The electronic band gap cannot be defined negative.')
+            return None
+
+        if not isinstance(self.value.magnitude, np.ndarray):  # for scalars
+            value = value[0]
+        return value * self.value.u
+
+    def resolve_type(self, logger: BoundLogger) -> Optional[str]:
+        """
+        Resolves the `type` of the electronic band gap based on the stored `momentum_transfer` values.
+
+        Args:
+            logger (BoundLogger): The logger to log messages.
+
+        Returns:
+            (Optional[str]): The resolved `type` of the electronic band gap.
+        """
+        mtr = self.momentum_transfer if self.momentum_transfer is not None else []
+
+        # Check if the `momentum_transfer` is [], and return the type and a warning in the log for `indirect` band gaps
+        if len(mtr) == 0:
+            if self.type == 'indirect':
+                logger.warning(
+                    'The `momentum_transfer` is not stored for an `indirect` band gap.'
+                )
+            return self.type
+
+        # Check if the `momentum_transfer` has at least two elements, and return None if it does not
+        if len(mtr) == 1:
+            logger.warning(
+                'The `momentum_transfer` should have at least two elements so that the difference can be calculated and the type of electronic band gap can be resolved.'
+            )
+            return None
+
+        # Resolve `type` from the difference between the initial and final momentum transfer
+        momentum_difference = np.diff(mtr, axis=0)
+        if (np.isclose(momentum_difference, np.zeros(3))).all():
+            return 'direct'
+        else:
+            return 'indirect'
+
+    def normalize(self, archive, logger) -> None:
+        super().normalize(archive, logger)
+
+        # Checks if the `value` is negative and sets it to None if it is.
+        self.value = self.validate_values(logger)
+        if self.value is None:
+            # ? What about deleting the class if `value` is None?
+            logger.error('The `value` of the electronic band gap is not stored.')
+            return
+
+        # Resolve the `type` of the electronic band gap from `momentum_transfer`, ONLY for scalar `value`
+        if isinstance(self.value.magnitude, np.ndarray):
+            logger.info(
+                'We do not support `type` which describe individual elements in an array `value`.'
+            )
+        else:
+            self.type = self.resolve_type(logger)
+
diff --git a/properties/band_structure.py b/properties/band_structure.py
new file mode 100644
index 00000000..239503be
--- /dev/null
+++ b/properties/band_structure.py
@@ -0,0 +1,328 @@
+#
+# Copyright The NOMAD Authors.
+#
+# This file is part of NOMAD. See https://nomad-lab.eu for further info.
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+#
+
+import numpy as np
+from structlog.stdlib import BoundLogger
+from typing import Optional, Tuple, Union
+import pint
+
+from nomad.metainfo import Quantity, Section, Context, SubSection
+
+from nomad_simulations.numerical_settings import KSpace
+from nomad_simulations.physical_property import (
+    PhysicalProperty,
+    validate_quantity_wrt_value,
+)
+from nomad_simulations.properties import ElectronicBandGap, FermiSurface
+from nomad_simulations.utils import get_sibling_section
+
+
+class BaseElectronicEigenvalues(PhysicalProperty):
+    """
+    A base section used to define basic quantities for the `ElectronicEigenvalues`  and `ElectronicBandStructure` properties.
+    """
+
+    iri = ''
+
+    n_bands = Quantity(
+        type=np.int32,
+        description="""
+        Number of bands / eigenvalues.
+        """,
+    )
+
+    value = Quantity(
+        type=np.float64,
+        unit='joule',
+        description="""
+        Value of the electronic eigenvalues.
+        """,
+    )
+
+    def __init__(
+        self, m_def: Section = None, m_context: Context = None, **kwargs
+    ) -> None:
+        super().__init__(m_def, m_context, **kwargs)
+        # ! `n_bands` need to be set up during initialization of the class
+        self.rank = [int(kwargs.get('n_bands'))]
+
+    def normalize(self, archive, logger) -> None:
+        super().normalize(archive, logger)
+
+
+class ElectronicEigenvalues(BaseElectronicEigenvalues):
+    """ """
+
+    iri = 'http://fairmat-nfdi.eu/taxonomy/ElectronicEigenvalues'
+
+    spin_channel = Quantity(
+        type=np.int32,
+        description="""
+        Spin channel of the corresponding electronic eigenvalues. It can take values of 0 or 1.
+        """,
+    )
+
+    occupation = Quantity(
+        type=np.float64,
+        shape=['*', 'n_bands'],
+        description="""
+        Occupation of the electronic eigenvalues. This is a number depending whether the `spin_channel` has been set or not.
+        If `spin_channel` is set, then this number is between 0 and 2, where 0 means that the state is unoccupied and 2 means
+        that the state is fully occupied; if `spin_channel` is not set, then this number is between 0 and 1. The shape of
+        this quantity is defined as `[K.n_points, K.dimensionality, n_bands]`, where `K` is a `variable` which can
+        be `KMesh` or `KLinePath`, depending whether the simulation mapped the whole Brillouin zone or just a specific
+        path.
+        """,
+    )
+
+    highest_occupied = Quantity(
+        type=np.float64,
+        unit='joule',
+        description="""
+        Highest occupied electronic eigenvalue. Together with `lowest_unoccupied`, it defines the
+        electronic band gap.
+        """,
+    )
+
+    lowest_unoccupied = Quantity(
+        type=np.float64,
+        unit='joule',
+        description="""
+        Lowest unoccupied electronic eigenvalue. Together with `highest_occupied`, it defines the
+        electronic band gap.
+        """,
+    )
+
+    # ? Should we add functionalities to handle min/max of the `value` in some specific cases, .e.g, bands around the Fermi level,
+    # ? core bands separated by gaps, and equivalently, higher-energy valence bands separated by gaps?
+
+    value_contributions = SubSection(
+        sub_section=BaseElectronicEigenvalues.m_def,
+        repeats=True,
+        description="""
+        Contributions to the electronic eigenvalues. Example, in the case of a DFT+GW calculation, the GW eigenvalues
+        are stored under `value`, and each contribution is identified by `label`:
+            - `'KS'`: Kohn-Sham contribution. This is also stored in the DFT entry under `ElectronicEigenvalues.value`.
+            - `'KSxc'`: Diagonal matrix elements of the expectation value of the Kohn-Sahm exchange-correlation potential.
+            - `'SigX'`: Diagonal matrix elements of the exchange self-energy. This is also stored in the GW entry under `ElectronicSelfEnergy.value`.
+            - `'SigC'`: Diagonal matrix elements of the correlation self-energy. This is also stored in the GW entry under `ElectronicSelfEnergy.value`.
+            - `'Zk'`: Quasiparticle renormalization factors contribution. This is also stored in the GW entry under `QuasiparticleWeights.value`.
+        """,
+    )
+
+    reciprocal_cell = Quantity(
+        type=KSpace.reciprocal_lattice_vectors,
+        description="""
+        Reference to the reciprocal lattice vectors stored under `KSpace`.
+        """,
+    )
+
+    def __init__(
+        self, m_def: Section = None, m_context: Context = None, **kwargs
+    ) -> None:
+        super().__init__(m_def, m_context, **kwargs)
+        self.name = self.m_def.name
+
+    @validate_quantity_wrt_value(name='occupation')
+    def order_eigenvalues(self) -> Union[bool, Tuple[pint.Quantity, np.ndarray]]:
+        """
+        Order the eigenvalues based on the `value` and `occupation`. The return `value` and
+        `occupation` are flattened.
+
+        Returns:
+            (Union[bool, Tuple[pint.Quantity, np.ndarray]]): The flattened and sorted `value` and `occupation`. If validation
+            fails, then it returns `False`.
+        """
+        total_shape = np.prod(self.value.shape)
+
+        # Order the indices in the flattened list of `value`
+        flattened_value = self.value.reshape(total_shape)
+        flattened_occupation = self.occupation.reshape(total_shape)
+        sorted_indices = np.argsort(flattened_value, axis=0)
+
+        sorted_value = (
+            np.take_along_axis(flattened_value.magnitude, sorted_indices, axis=0)
+            * flattened_value.u
+        )
+        sorted_occupation = np.take_along_axis(
+            flattened_occupation, sorted_indices, axis=0
+        )
+        self.m_cache['sorted_eigenvalues'] = True
+        return sorted_value, sorted_occupation
+
+    def resolve_homo_lumo_eigenvalues(
+        self,
+    ) -> Tuple[Optional[pint.Quantity], Optional[pint.Quantity]]:
+        """
+        Resolve the `highest_occupied` and `lowest_unoccupied` eigenvalues by performing a binary search on the
+        flattened and sorted `value` and `occupation`. If these quantities already exist, overwrite them or return
+        them if it is not possible to resolve from `value` and `occupation`.
+
+        Returns:
+            (Tuple[Optional[pint.Quantity], Optional[pint.Quantity]]): The `highest_occupied` and
+            `lowest_unoccupied` eigenvalues.
+        """
+        # Sorting `value` and `occupation`
+        if not self.order_eigenvalues():  # validation fails
+            if self.highest_occupied is not None and self.lowest_unoccupied is not None:
+                return self.highest_occupied, self.lowest_unoccupied
+            return None, None
+        sorted_value, sorted_occupation = self.order_eigenvalues()
+        sorted_value_unit = sorted_value.u
+        sorted_value = sorted_value.magnitude
+
+        # Binary search ot find the transition point between `occupation = 2` and `occupation = 0`
+        tolerance = 1e-3  # TODO add tolerance from config fields
+        homo = self.highest_occupied
+        lumo = self.lowest_unoccupied
+        mid = np.searchsorted(sorted_occupation <= tolerance, True) - 1
+        if mid >= 0 and mid < len(sorted_occupation) - 1:
+            if sorted_occupation[mid] > 0 and (
+                sorted_occupation[mid + 1] >= -tolerance
+                and sorted_occupation[mid + 1] <= tolerance
+            ):
+                homo = sorted_value[mid] * sorted_value_unit
+                lumo = sorted_value[mid + 1] * sorted_value_unit
+
+        return homo, lumo
+
+    def extract_band_gap(self) -> Optional[ElectronicBandGap]:
+        """
+        Extract the electronic band gap from the `highest_occupied` and `lowest_unoccupied` eigenvalues.
+        If the difference of `highest_occupied` and `lowest_unoccupied` is negative, the band gap `value` is set to 0.0.
+
+        Returns:
+            (Optional[ElectronicBandGap]): The extracted electronic band gap section to be stored in `Outputs`.
+        """
+        band_gap = None
+        homo, lumo = self.resolve_homo_lumo_eigenvalues()
+        if homo and lumo:
+            band_gap = ElectronicBandGap(is_derived=True, physical_property_ref=self)
+
+            if (lumo - homo).magnitude < 0:
+                band_gap.value = 0.0
+            else:
+                band_gap.value = lumo - homo
+        return band_gap
+
+    # TODO fix this method once `FermiSurface` property is implemented
+    def extract_fermi_surface(self, logger: BoundLogger) -> Optional[FermiSurface]:
+        """
+        Extract the Fermi surface for metal systems and using the `FermiLevel.value`.
+
+        Args:
+            logger (BoundLogger): The logger to log messages.
+
+        Returns:
+            (Optional[FermiSurface]): The extracted Fermi surface section to be stored in `Outputs`.
+        """
+        # Check if the system has a finite band gap
+        homo, lumo = self.resolve_homo_lumo_eigenvalues()
+        if (homo and lumo) and (lumo - homo).magnitude > 0:
+            return None
+
+        # Get the `fermi_level.value`
+        fermi_level = get_sibling_section(
+            section=self, sibling_section_name='fermi_level', logger=logger
+        )
+        if fermi_level is None:
+            logger.warning(
+                'Could not extract the `FermiSurface`, because `FermiLevel` is not stored.'
+            )
+            return None
+        fermi_level_value = fermi_level.value.magnitude
+
+        # Extract values close to the `fermi_level.value`
+        tolerance = 1e-8  # TODO change this for a config field
+        fermi_indices = np.logical_and(
+            self.value.magnitude >= (fermi_level_value - tolerance),
+            self.value.magnitude <= (fermi_level_value + tolerance),
+        )
+        fermi_values = self.value[fermi_indices]
+
+        # Store `FermiSurface` values
+        # ! This is wrong (!) the `value` should be the `KMesh.points`, not the `ElectronicEigenvalues.value`
+        fermi_surface = FermiSurface(
+            n_bands=self.n_bands,
+            is_derived=True,
+            physical_property_ref=self,
+        )
+        fermi_surface.variables = self.variables
+        fermi_surface.value = fermi_values
+        return fermi_surface
+
+    def resolve_reciprocal_cell(self) -> Optional[pint.Quantity]:
+        """
+        Resolve the reciprocal cell from the `KSpace` numerical settings section.
+
+        Returns:
+            Optional[pint.Quantity]: _description_
+        """
+        numerical_settings = self.m_xpath(
+            'm_parent.m_parent.model_method[-1].numerical_settings', dict=False
+        )
+        if numerical_settings is None:
+            return None
+        k_space = None
+        for setting in numerical_settings:
+            if isinstance(setting, KSpace):
+                k_space = setting
+                break
+        if k_space is None:
+            return None
+        return k_space
+
+    def normalize(self, archive, logger) -> None:
+        super().normalize(archive, logger)
+
+        # Resolve `highest_occupied` and `lowest_unoccupied` eigenvalues
+        self.highest_occupied, self.lowest_unoccupied = (
+            self.resolve_homo_lumo_eigenvalues()
+        )
+
+        # `ElectronicBandGap` extraction
+        band_gap = self.extract_band_gap()
+        if band_gap is not None:
+            self.m_parent.electronic_band_gaps.append(band_gap)
+
+        # TODO uncomment once `FermiSurface` property is implemented
+        # `FermiSurface` extraction
+        # fermi_surface = self.extract_fermi_surface(logger)
+        # if fermi_surface is not None:
+        #     self.m_parent.fermi_surfaces.append(fermi_surface)
+
+        # Resolve `reciprocal_cell` from the `KSpace` numerical settings section
+        self.reciprocal_cell = self.resolve_reciprocal_cell()
+
+
+class ElectronicBandStructure(ElectronicEigenvalues):
+    """
+    Accessible energies by the charges (electrons and holes) in the reciprocal space.
+    """
+
+    iri = 'http://fairmat-nfdi.eu/taxonomy/ElectronicBandStructure'
+
+    def __init__(
+        self, m_def: Section = None, m_context: Context = None, **kwargs
+    ) -> None:
+        super().__init__(m_def, m_context, **kwargs)
+        self.name = self.m_def.name
+
+    def normalize(self, archive, logger) -> None:
+        super().normalize(archive, logger)
+
diff --git a/properties/energies.py b/properties/energies.py
index 800f69d4..8ef7c7fb 100644
--- a/properties/energies.py
+++ b/properties/energies.py
@@ -38,6 +38,85 @@
 from nomad.metainfo import Quantity, SubSection, MEnum
 from .physical_property import PhysicalProperty
 
+#
+# Copyright The NOMAD Authors.
+#
+# This file is part of NOMAD. See https://nomad-lab.eu for further info.
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+#
+
+import numpy as np
+
+from nomad.metainfo import Quantity, Section, Context
+
+from nomad_simulations.physical_property import PhysicalProperty
+
+
+class FermiLevel(PhysicalProperty):
+    """
+    Energy required to add or extract a charge from a material at zero temperature. It can be also defined as the chemical potential at zero temperature.
+    """
+
+    # ! implement `iri` and `rank` as part of `m_def = Section()`
+
+    iri = 'http://fairmat-nfdi.eu/taxonomy/FermiLevel'
+
+    value = Quantity(
+        type=np.float64,
+        unit='joule',
+        description="""
+        The value of the Fermi level.
+        """,
+    )
+
+    def __init__(
+        self, m_def: Section = None, m_context: Context = None, **kwargs
+    ) -> None:
+        super().__init__(m_def, m_context, **kwargs)
+        self.rank = []
+        self.name = self.m_def.name
+
+    def normalize(self, archive, logger) -> None:
+        super().normalize(archive, logger)
+
+
+class ChemicalPotential(PhysicalProperty):
+    """
+    Free energy cost of adding or extracting a particle from a thermodynamic system.
+    """
+
+    # ! implement `iri` and `rank` as part of `m_def = Section()`
+
+    iri = 'http://fairmat-nfdi.eu/taxonomy/ChemicalPotential'
+
+    value = Quantity(
+        type=np.float64,
+        unit='joule',
+        description="""
+        The value of the chemical potential.
+        """,
+    )
+
+    def __init__(
+        self, m_def: Section = None, m_context: Context = None, **kwargs
+    ) -> None:
+        super().__init__(m_def, m_context, **kwargs)
+        self.rank = []
+        self.name = self.m_def.name
+
+    def normalize(self, archive, logger) -> None:
+        super().normalize(archive, logger)
 
 
 class EnergyContributions(PhysicalProperty):
diff --git a/properties/fermi_surface.py b/properties/fermi_surface.py
new file mode 100644
index 00000000..4b473a9a
--- /dev/null
+++ b/properties/fermi_surface.py
@@ -0,0 +1,53 @@
+#
+# Copyright The NOMAD Authors.
+#
+# This file is part of NOMAD. See https://nomad-lab.eu for further info.
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+#
+
+import numpy as np
+
+from nomad.metainfo import Quantity, Section, Context
+
+from nomad_simulations.physical_property import PhysicalProperty
+
+
+# TODO This class is not implemented yet. @JosePizarro3 will work in another PR to implement it.
+class FermiSurface(PhysicalProperty):
+    """
+    Energy boundary in reciprocal space that separates the filled and empty electronic states in a metal.
+    It is related with the crossing points in reciprocal space by the chemical potential or, equivalently at
+    zero temperature, the Fermi level.
+    """
+
+    iri = 'http://fairmat-nfdi.eu/taxonomy/FermiSurface'
+
+    n_bands = Quantity(
+        type=np.int32,
+        description="""
+        Number of bands / eigenvalues.
+        """,
+    )
+
+    def __init__(
+        self, m_def: Section = None, m_context: Context = None, **kwargs
+    ) -> None:
+        super().__init__(m_def, m_context, **kwargs)
+        # ! `n_bands` need to be set up during initialization of the class
+        self.rank = [int(kwargs.get('n_bands'))]
+        self.name = self.m_def.name
+
+    def normalize(self, archive, logger) -> None:
+        super().normalize(archive, logger)
+
diff --git a/properties/hopping_matrix.py b/properties/hopping_matrix.py
new file mode 100644
index 00000000..cebeea71
--- /dev/null
+++ b/properties/hopping_matrix.py
@@ -0,0 +1,105 @@
+#
+# Copyright The NOMAD Authors.
+#
+# This file is part of NOMAD. See https://nomad-lab.eu for further info.
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+#
+
+import numpy as np
+
+from nomad.metainfo import Quantity, Section, Context
+
+from nomad_simulations.physical_property import PhysicalProperty
+
+
+class HoppingMatrix(PhysicalProperty):
+    """
+    Transition probability between two atomic orbitals in a tight-binding model.
+    """
+
+    iri = 'http://fairmat-nfdi.eu/taxonomy/HoppingMatrix'
+
+    n_orbitals = Quantity(
+        type=np.int32,
+        description="""
+        Number of orbitals in the tight-binding model.
+        """,
+    )
+
+    degeneracy_factors = Quantity(
+        type=np.int32,
+        shape=['*'],
+        description="""
+        Degeneracy of each Wigner-Seitz point.
+        """,
+    )
+
+    value = Quantity(
+        type=np.complex128,
+        unit='joule',
+        description="""
+        Value of the hopping matrix in joules. The elements are complex numbers defined for each Wigner-Seitz point and
+        each pair of orbitals; thus, `rank = [n_orbitals, n_orbitals]`. Note this contains also the onsite values, i.e.,
+        it includes the Wigner-Seitz point (0, 0, 0), hence the `CrystalFieldSplitting` values.
+        """,
+    )
+
+    def __init__(
+        self, m_def: Section = None, m_context: Context = None, **kwargs
+    ) -> None:
+        super().__init__(m_def, m_context, **kwargs)
+        # ! n_orbitals need to be set up during initialization of the class
+        self.rank = [int(kwargs.get('n_orbitals')), int(kwargs.get('n_orbitals'))]
+        self.name = self.m_def.name
+
+    # TODO add normalization to extract DOS, band structure, etc, properties from `HoppingMatrix`
+
+    def normalize(self, archive, logger) -> None:
+        super().normalize(archive, logger)
+
+
+class CrystalFieldSplitting(PhysicalProperty):
+    """
+    Energy difference between the degenerated orbitals of an ion in a crystal field environment.
+    """
+
+    iri = 'http://fairmat-nfdi.eu/taxonomy/CrystalFieldSplitting'
+
+    n_orbitals = Quantity(
+        type=np.int32,
+        description="""
+        Number of orbitals in the tight-binding model.
+        """,
+    )
+
+    value = Quantity(
+        type=np.float64,
+        unit='joule',
+        description="""
+        Value of the crystal field splittings in joules. This is the intra-orbital local contribution, i.e., the same orbital
+        at the same Wigner-Seitz point (0, 0, 0).
+        """,
+    )
+
+    def __init__(
+        self, m_def: Section = None, m_context: Context = None, **kwargs
+    ) -> None:
+        super().__init__(m_def, m_context, **kwargs)
+        # ! `n_orbitals` need to be set up during initialization of the class
+        self.rank = [int(kwargs.get('n_orbitals'))]
+        self.name = self.m_def.name
+
+    def normalize(self, archive, logger) -> None:
+        super().normalize(archive, logger)
+
diff --git a/properties/permittivity.py b/properties/permittivity.py
new file mode 100644
index 00000000..d97c5cd4
--- /dev/null
+++ b/properties/permittivity.py
@@ -0,0 +1,119 @@
+#
+# Copyright The NOMAD Authors.
+#
+# This file is part of NOMAD. See https://nomad-lab.eu for further info.
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+#
+
+import numpy as np
+from structlog.stdlib import BoundLogger
+from typing import Optional, List
+
+from nomad.metainfo import Quantity, Section, Context, MEnum
+
+from nomad_simulations.physical_property import PhysicalProperty
+from nomad_simulations.utils import get_variables
+from nomad_simulations.variables import Frequency, KMesh
+from nomad_simulations.properties.spectral_profile import AbsorptionSpectrum
+
+
+# TODO add `DielectricStrength` when we have examples and understand how to extract it from the `Permittivity` tensor.
+
+
+class Permittivity(PhysicalProperty):
+    """
+    Response of the material to polarize in the presence of an electric field.
+
+    Alternative names: `DielectricFunction`.
+    """
+
+    iri = 'http://fairmat-nfdi.eu/taxonomy/Permittivity'
+
+    type = Quantity(
+        type=MEnum('static', 'dynamic'),
+        description="""
+        Type of permittivity which allows to identify if the permittivity depends on the frequency or not.
+        """,
+    )
+
+    value = Quantity(
+        type=np.complex128,
+        # unit='joule',  # TODO check units (they have to match `SpectralProfile.value`)
+        description="""
+        Value of the permittivity tensor. If the value does not depend on the scattering vector `q`, then we
+        can extract the optical absorption spectrum from the imaginary part of the permittivity tensor (this is also called
+        macroscopic dielectric function).
+        """,
+    )
+
+    # ? We need use cases to understand if we need to define contributions to the permittivity tensor.
+    # ? `ionic` and `electronic` contributions are common in the literature.
+
+    def __init__(
+        self, m_def: Section = None, m_context: Context = None, **kwargs
+    ) -> None:
+        super().__init__(m_def, m_context, **kwargs)
+        self.rank = [3, 3]
+        self.name = self.m_def.name
+        self._axes_map = ['xx', 'yy', 'zz']
+
+    def resolve_type(self) -> str:
+        frequencies = get_variables(self.variables, Frequency)
+        if len(frequencies) > 0:
+            return 'dynamic'
+        return 'static'
+
+    def extract_absorption_spectra(
+        self, logger: BoundLogger
+    ) -> Optional[List[AbsorptionSpectrum]]:
+        """
+        Extract the absorption spectrum from the imaginary part of the permittivity tensor.
+        """
+        # If the `pemittivity` depends on the scattering vector `q`, then we cannot extract the absorption spectrum
+        q_mesh = get_variables(self.variables, KMesh)
+        if len(q_mesh) > 0:
+            logger.warning(
+                'The `permittivity` depends on the scattering vector `q`, so that we cannot extract the absorption spectrum.'
+            )
+            return None
+        # Extract the `Frequency` variable to extract the absorption spectrum
+        frequencies = get_variables(self.variables, Frequency)
+        if len(frequencies) == 0:
+            logger.warning(
+                'The `permittivity` does not have a `Frequency` variable to extract the absorption spectrum.'
+            )
+            return None
+        # Define the `absorption_spectra` for each principal direction along the diagonal of the `Permittivity.value` as the imaginary part
+        spectra = []
+        for i in range(3):
+            val = self.value[:, i, i].imag
+            absorption_spectrum = AbsorptionSpectrum(
+                axis=self._axes_map[i], variables=frequencies
+            )
+            absorption_spectrum.value = val
+            absorption_spectrum.physical_property_ref = self
+            spectra.append(absorption_spectrum)
+        return spectra
+
+    def normalize(self, archive, logger) -> None:
+        super().normalize(archive, logger)
+
+        # Resolve the `type` of permittivity
+        self.type = self.resolve_type()
+
+        # `AbsorptionSpectrum` extraction
+        absorption_spectra = self.extract_absorption_spectra(logger)
+        if absorption_spectra is not None:
+            self.m_parent.absorption_spectrum = absorption_spectra
+
diff --git a/properties/spectral_profile.py b/properties/spectral_profile.py
new file mode 100644
index 00000000..f2c894f0
--- /dev/null
+++ b/properties/spectral_profile.py
@@ -0,0 +1,619 @@
+#
+# Copyright The NOMAD Authors.
+#
+# This file is part of NOMAD. See https://nomad-lab.eu for further info.
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+#
+
+import numpy as np
+from structlog.stdlib import BoundLogger
+from typing import Optional, List, Dict
+import pint
+
+from nomad import config
+from nomad.metainfo import Quantity, SubSection, Section, Context, MEnum
+
+from nomad_simulations.utils import get_sibling_section, get_variables
+from nomad_simulations.physical_property import PhysicalProperty
+from nomad_simulations.variables import Energy2 as Energy
+from nomad_simulations.atoms_state import AtomsState, OrbitalsState
+from nomad_simulations.properties.band_gap import ElectronicBandGap
+
+
+class SpectralProfile(PhysicalProperty):
+    """
+    A base section used to define the spectral profile.
+    """
+
+    value = Quantity(
+        type=np.float64,
+        description="""
+        The value of the intensities of a spectral profile in arbitrary units.
+        """,
+    )  # TODO check units and normalization_factor of DOS and Spectras and see whether they can be merged
+
+    def __init__(
+        self, m_def: Section = None, m_context: Context = None, **kwargs
+    ) -> None:
+        super().__init__(m_def, m_context, **kwargs)
+        self.rank = []
+
+    def is_valid_spectral_profile(self) -> bool:
+        """
+        Check if the spectral profile is valid, i.e., if all `value` are defined positive.
+
+        Returns:
+            (bool): True if the spectral profile is valid, False otherwise.
+        """
+        if (self.value < 0.0).any():
+            return False
+        return True
+
+    def normalize(self, archive, logger) -> None:
+        super().normalize(archive, logger)
+
+        if self.is_valid_spectral_profile() is False:
+            logger.error(
+                'Invalid negative intensities found: could not validate spectral profile.'
+            )
+            return
+
+
+class DOSProfile(SpectralProfile):
+    """
+    A base section used to define the `value` of the `ElectronicDensityOfState` property. This is useful when containing
+    contributions for `projected_dos` with the correct unit.
+    """
+
+    value = Quantity(
+        type=np.float64,
+        unit='1/joule',
+        description="""
+        The value of the electronic DOS.
+        """,
+    )
+
+    def __init__(
+        self, m_def: Section = None, m_context: Context = None, **kwargs
+    ) -> None:
+        super().__init__(m_def, m_context, **kwargs)
+
+    def resolve_pdos_name(self, logger: BoundLogger) -> Optional[str]:
+        """
+        Resolve the `name` of the projected `DOSProfile` from the `entity_ref` section. This is resolved as:
+            - `'atom X'` with 'X' being the chemical symbol for `AtomsState` references.
+            -  `'orbital Y X'` with 'X' being the chemical symbol and 'Y' the orbital label for `OrbitalsState` references.
+
+        Args:
+            logger (BoundLogger): The logger to log messages.
+
+        Returns:
+            (Optional[str]): The resolved `name` of the projected DOS profile.
+        """
+        if self.entity_ref is None:
+            logger.warning(
+                'The `entity_ref` is not set for the DOS profile. Could not resolve the `name`.'
+            )
+            return None
+
+        # Resolve the `name` from the `entity_ref`
+        name = None
+        if isinstance(self.entity_ref, AtomsState):
+            name = f'atom {self.entity_ref.chemical_symbol}'
+        elif isinstance(self.entity_ref, OrbitalsState):
+            name = f'orbital {self.entity_ref.l_quantum_symbol}{self.entity_ref.ml_quantum_symbol} {self.entity_ref.m_parent.chemical_symbol}'
+        return name
+
+    def normalize(self, archive, logger) -> None:
+        super().normalize(archive, logger)
+
+        # We resolve
+        self.name = self.resolve_pdos_name(logger)
+
+
+class ElectronicDensityOfStates(DOSProfile):
+    """
+    Number of electronic states accessible for the charges per energy and per volume.
+    """
+
+    # ! implement `iri` and `rank` as part of `m_def = Section()`
+    iri = 'http://fairmat-nfdi.eu/taxonomy/ElectronicDensityOfStates'
+
+    spin_channel = Quantity(
+        type=np.int32,
+        description="""
+        Spin channel of the corresponding electronic DOS. It can take values of 0 or 1.
+        """,
+    )
+
+    # TODO clarify the role of `energies_origin` once `ElectronicEigenvalues` is implemented
+    energies_origin = Quantity(
+        type=np.float64,
+        unit='joule',
+        description="""
+        Energy level denoting the origin along the energy axis, used for comparison and visualization. It is
+        defined as the `ElectronicEigenvalues.highest_occupied_energy`.
+        """,
+    )
+
+    normalization_factor = Quantity(
+        type=np.float64,
+        description="""
+        Normalization factor for electronic DOS to get a cell-independent intensive DOS. The cell-independent
+        intensive DOS is as the integral from the lowest (most negative) energy to the Fermi level for a neutrally
+        charged system (i.e., the sum of `AtomsState.charge` is zero).
+        """,
+    )
+
+    # ? Do we want to store the integrated value here os as part of an nomad-analysis tool? Check `dos_integrator.py` module in dos normalizer repository
+    # value_integrated = Quantity(
+    #     type=np.float64,
+    #     description="""
+    #     The cumulative intensities integrated from from the lowest (most negative) energy to the Fermi level.
+    #     """,
+    # )
+
+    projected_dos = SubSection(
+        sub_section=DOSProfile.m_def,
+        repeats=True,
+        description="""
+        Projected DOS. It can be atom- (different elements in the unit cell) or orbital-projected. These can be calculated in a cascade as:
+            - If the total DOS is not present, we sum all atom-projected DOS to obtain it.
+            - If the atom-projected DOS is not present, we sum all orbital-projected DOS to obtain it.
+        Note: the cover given by summing up contributions is not perfect, and will depend on the projection functions used.
+
+        In `projected_dos`, `name` and `entity_ref` must be set in order for normalization to work:
+            - The `entity_ref` is the `OrbitalsState` or `AtomsState` sections.
+            - The `name` of the projected DOS should be `'atom X'` or `'orbital Y X'`, with 'X' being the chemical symbol and 'Y' the orbital label.
+            These can be extracted from `entity_ref`.
+        """,
+    )
+
+    def __init__(
+        self, m_def: Section = None, m_context: Context = None, **kwargs
+    ) -> None:
+        super().__init__(m_def, m_context, **kwargs)
+        self.name = self.m_def.name
+
+    def resolve_energies_origin(
+        self,
+        energies: pint.Quantity,
+        fermi_level: Optional[pint.Quantity],
+        logger: BoundLogger,
+    ) -> Optional[pint.Quantity]:
+        """
+        Resolve the origin of reference for the energies from the sibling `ElectronicEigenvalues` section and its
+        `highest_occupied` level, or if this does not exist, from the `fermi_level` value as extracted from the sibling property, `FermiLevel`.
+
+        Args:
+            fermi_level (Optional[pint.Quantity]): The resolved Fermi level.
+            energies (pint.Quantity): The grid points of the `Energy` variable.
+            logger (BoundLogger): The logger to log messages.
+
+        Returns:
+            (Optional[pint.Quantity]): The resolved origin of reference for the energies.
+        """
+        # Check if the variables contain more than one variable (different than Energy)
+        # ? Is this correct or should be use the index of energies to extract the proper shape element in `self.value` being used for `dos_values`?
+        if len(self.variables) > 1:
+            logger.warning(
+                'The ElectronicDensityOfStates section contains more than one variable. We cannot extract the energy reference.'
+            )
+            return None
+
+        # Extract the `ElectronicEigenvalues` section to get the `highest_occupied` and `lowest_unoccupied` energies
+        # TODO implement once `ElectronicEigenvalues` is in the schema
+        eigenvalues = get_sibling_section(
+            section=self, sibling_section_name='electronic_eigenvalues', logger=logger
+        )  # we consider `index_sibling` to be 0
+        highest_occupied_energy = (
+            eigenvalues.highest_occupied if eigenvalues is not None else None
+        )
+        lowest_unoccupied_energy = (
+            eigenvalues.lowest_unoccupied if eigenvalues is not None else None
+        )
+        # and set defaults for `highest_occupied_energy` and `lowest_unoccupied_energy` in `m_cache`
+        if highest_occupied_energy is not None:
+            self.m_cache['highest_occupied_energy'] = highest_occupied_energy
+        if lowest_unoccupied_energy is not None:
+            self.m_cache['lowest_unoccupied_energy'] = lowest_unoccupied_energy
+
+        # Set thresholds for the energies and values
+        energy_threshold = config.normalize.band_structure_energy_tolerance
+        value_threshold = 1e-8  # The DOS value that is considered to be zero
+
+        # Check that the closest `energies` to the energy reference is not too far away.
+        # If it is very far away, normalization may be very inaccurate and we do not report it.
+        dos_values = self.value.magnitude
+        eref = highest_occupied_energy if fermi_level is None else fermi_level
+        fermi_idx = (np.abs(energies - eref)).argmin()
+        fermi_energy_closest = energies[fermi_idx]
+        distance = np.abs(fermi_energy_closest - eref)
+        single_peak_fermi = False
+        if distance.magnitude <= energy_threshold:
+            # See if there are zero values close below the energy reference.
+            idx = fermi_idx
+            idx_descend = fermi_idx
+            while True:
+                try:
+                    value = dos_values[idx]
+                    energy_distance = np.abs(eref - energies[idx])
+                except IndexError:
+                    break
+                if energy_distance.magnitude > energy_threshold:
+                    break
+                if value <= value_threshold:
+                    idx_descend = idx
+                    break
+                idx -= 1
+
+            # See if there are zero values close above the fermi energy.
+            idx = fermi_idx
+            idx_ascend = fermi_idx
+            while True:
+                try:
+                    value = dos_values[idx]
+                    energy_distance = np.abs(eref - energies[idx])
+                except IndexError:
+                    break
+                if energy_distance.magnitude > energy_threshold:
+                    break
+                if value <= value_threshold:
+                    idx_ascend = idx
+                    break
+                idx += 1
+
+            # If there is a single peak at fermi energy, no
+            # search needs to be performed.
+            if idx_ascend != fermi_idx and idx_descend != fermi_idx:
+                self.m_cache['highest_occupied_energy'] = fermi_energy_closest
+                self.m_cache['lowest_unoccupied_energy'] = fermi_energy_closest
+                single_peak_fermi = True
+
+            if not single_peak_fermi:
+                # Look for highest occupied energy below the descend index
+                idx = idx_descend
+                while True:
+                    try:
+                        value = dos_values[idx]
+                    except IndexError:
+                        break
+                    if value > value_threshold:
+                        idx = idx if idx == idx_descend else idx + 1
+                        self.m_cache['highest_occupied_energy'] = energies[idx]
+                        break
+                    idx -= 1
+                # Look for lowest unoccupied energy above idx_ascend
+                idx = idx_ascend
+                while True:
+                    try:
+                        value = dos_values[idx]
+                    except IndexError:
+                        break
+                    if value > value_threshold:
+                        idx = idx if idx == idx_ascend else idx - 1
+                        self.m_cache['highest_occupied_energy'] = energies[idx]
+                        break
+                    idx += 1
+
+        # Return the `highest_occupied_energy` as the `energies_origin`, or the `fermi_level` if it is not None
+        energies_origin = self.m_cache.get('highest_occupied_energy')
+        if energies_origin is None:
+            energies_origin = fermi_level
+        return energies_origin
+
+    def resolve_normalization_factor(self, logger: BoundLogger) -> Optional[float]:
+        """
+        Resolve the `normalization_factor` for the electronic DOS to get a cell-independent intensive DOS.
+
+        Args:
+            logger (BoundLogger): The logger to log messages.
+
+        Returns:
+            (Optional[float]): The normalization factor.
+        """
+        # Get the `ModelSystem` as referenced in the `Outputs.model_system_ref`
+        model_system = get_sibling_section(
+            section=self, sibling_section_name='model_system_ref', logger=logger
+        )
+        if model_system is None:
+            logger.warning(
+                'Could not resolve the referenced `ModelSystem` in the `Outputs`.'
+            )
+            return None
+
+        # Get the originally parsed `AtomicCell`, which is the first element stored in `ModelSystem.cell` of name `'AtomicCell'`
+        atomic_cell = None
+        for cell in model_system.cell:
+            if cell.name == 'AtomicCell':  # we get the originally parsed `AtomicCell`
+                atomic_cell = cell
+                break
+        if atomic_cell is None:
+            logger.warning(
+                'Could not resolve the `AtomicCell` from the referenced `ModelSystem`.'
+            )
+            return None
+
+        # Get the `atoms_state` and their `atomic_number` from the `AtomicCell`
+        if atomic_cell.atoms_state is None or len(atomic_cell.atoms_state) == 0:
+            logger.warning('Could not resolve the `atoms_state` from the `AtomicCell`.')
+            return None
+        atomic_numbers = [atom.atomic_number for atom in atomic_cell.atoms_state]
+
+        # Return `normalization_factor` depending if the calculation is spin polarized or not
+        if self.spin_channel is not None:
+            normalization_factor = 1 / (2 * sum(atomic_numbers))
+        else:
+            normalization_factor = 1 / sum(atomic_numbers)
+        return normalization_factor
+
+    def extract_band_gap(self) -> Optional[ElectronicBandGap]:
+        """
+        Extract the electronic band gap from the `highest_occupied_energy` and `lowest_unoccupied_energy` stored
+        in `m_cache` from `resolve_energies_origin()`. If the difference of `highest_occupied_energy` and
+        `lowest_unoccupied_energy` is negative, the band gap `value` is set to 0.0.
+
+        Returns:
+            (Optional[ElectronicBandGap]): The extracted electronic band gap section to be stored in `Outputs`.
+        """
+        band_gap = None
+        homo = self.m_cache.get('highest_occupied_energy')
+        lumo = self.m_cache.get('lowest_unoccupied_energy')
+        if homo and lumo:
+            band_gap = ElectronicBandGap()
+            band_gap.is_derived = True
+            band_gap.physical_property_ref = self
+
+            if (homo - lumo).magnitude < 0:
+                band_gap.value = 0.0
+            else:
+                band_gap.value = homo - lumo
+        return band_gap
+
+    def extract_projected_dos(
+        self, type: str, logger: BoundLogger
+    ) -> List[Optional[DOSProfile]]:
+        """
+        Extract the projected DOS from the `projected_dos` section and the specified `type`.
+
+        Args:
+            type (str): The type of the projected DOS to extract. It can be `'atom'` or `'orbital'`.
+
+        Returns:
+            (DOSProfile): The extracted projected DOS.
+        """
+        extracted_pdos = []
+        for pdos in self.projected_dos:
+            # We make sure each PDOS is normalized
+            pdos.normalize(None, logger)
+
+            # Initial check for `name` and `entity_ref`
+            if (
+                pdos.name is None
+                or pdos.entity_ref is None
+                or len(pdos.entity_ref) == 0
+            ):
+                logger.warning(
+                    '`name` or `entity_ref` are not set for `projected_dos` and they are required for normalization to work.'
+                )
+                return None
+
+            if type in pdos.name:
+                extracted_pdos.append(pdos)
+        return extracted_pdos
+
+    def generate_from_projected_dos(
+        self, logger: BoundLogger
+    ) -> Optional[pint.Quantity]:
+        """
+        Generate the total `value` of the electronic DOS from the `projected_dos` contributions. If the `projected_dos`
+        is not present, it returns `None`.
+
+        Args:
+            logger (BoundLogger): The logger to log messages.
+
+        Returns:
+            (Optional[pint.Quantity]): The total `value` of the electronic DOS.
+        """
+        if self.projected_dos is None or len(self.projected_dos) == 0:
+            return None
+
+        # Extract `Energy` variables
+        energies = get_variables(self.variables, Energy)
+        if len(energies) != 1:
+            logger.warning(
+                'The `ElectronicDensityOfStates` does not contain an `Energy` variable to extract the DOS.'
+            )
+            return None
+
+        # We distinguish between orbital and atom `projected_dos`
+        orbital_projected = self.extract_projected_dos('orbital', logger)
+        atom_projected = self.extract_projected_dos('atom', logger)
+
+        # Extract `atom_projected` from `orbital_projected` by summing up the `orbital_projected` contributions for each atom
+        if len(atom_projected) == 0:
+            atom_data: Dict[AtomsState, List[DOSProfile]] = {}
+            for orb_pdos in orbital_projected:
+                # `entity_ref` is the `OrbitalsState` section, whose parent is `AtomsState`
+                entity_ref = orb_pdos.entity_ref.m_parent
+                if entity_ref in atom_data:
+                    atom_data[entity_ref].append(orb_pdos)
+                else:
+                    atom_data[entity_ref] = [orb_pdos]
+            for ref, data in atom_data.items():
+                atom_dos = DOSProfile(
+                    name=f'atom {ref.chemical_symbol}',
+                    entity_ref=ref,
+                    variables=energies,
+                )
+                orbital_values = [
+                    dos.value.magnitude for dos in data
+                ]  # to avoid warnings from pint
+                orbital_unit = data[0].value.u
+                atom_dos.value = np.sum(orbital_values, axis=0) * orbital_unit
+                atom_projected.append(atom_dos)
+            # We concatenate the `atom_projected` to the `projected_dos`
+            self.projected_dos = orbital_projected + atom_projected
+
+        # Extract `value` from `atom_projected` by summing up the `atom_projected` contributions
+        value = self.value
+        if value is None:
+            atom_values = [
+                dos.value.magnitude for dos in atom_projected
+            ]  # to avoid warnings from pint
+            atom_unit = atom_projected[0].value.u
+            value = np.sum(atom_values, axis=0) * atom_unit
+        return value
+
+    def normalize(self, archive, logger) -> None:
+        super().normalize(archive, logger)
+
+        # Initial check to see if `variables` contains the required `Energy` variable
+        energies = get_variables(self.variables, Energy)
+        if len(energies) != 0:
+            return
+        energies = energies[0].points
+
+        # Resolve `fermi_level` from a sibling section with respect to `ElectronicDensityOfStates`
+        fermi_level = get_sibling_section(
+            section=self, sibling_section_name='fermi_level', logger=logger
+        )  # * we consider `index_sibling` to be 0
+        if fermi_level is not None:
+            fermi_level = fermi_level.value
+        # and the `energies_origin` from the sibling `ElectronicEigenvalues` section
+        self.energies_origin = self.resolve_energies_origin(
+            energies, fermi_level, logger
+        )
+        if self.energies_origin is None:
+            logger.info('Could not resolve the `energies_origin` for the DOS')
+
+        # Resolve `normalization_factor`
+        if self.normalization_factor is None:
+            self.normalization_factor = self.resolve_normalization_factor(logger)
+
+        # `ElectronicBandGap` extraction
+        band_gap = self.extract_band_gap()
+        if band_gap is not None:
+            self.m_parent.electronic_band_gap.append(band_gap)
+
+        # Total `value` extraction from `projected_dos`
+        value_from_pdos = self.generate_from_projected_dos(logger)
+        if self.value is None and value_from_pdos is not None:
+            logger.info(
+                'The `ElectronicDensityOfStates.value` is not stored. We will attempt to obtain it by summing up projected DOS contributions, if these are present.'
+            )
+            self.value = value_from_pdos
+
+
+class AbsorptionSpectrum(SpectralProfile):
+    """ """
+
+    # ! implement `iri` and `rank` as part of `m_def = Section()`
+
+    axis = Quantity(
+        type=MEnum('xx', 'yy', 'zz'),
+        description="""
+        Axis of the absorption spectrum. This is related with the polarization direction, and can be seen as the
+        principal term in the tensor `Permittivity.value` (see permittivity.py module).
+        """,
+    )
+
+    def __init__(
+        self, m_def: Section = None, m_context: Context = None, **kwargs
+    ) -> None:
+        super().__init__(m_def, m_context, **kwargs)
+        # Set the name of the section
+        self.name = self.m_def.name
+
+    def normalize(self, archive, logger) -> None:
+        super().normalize(archive, logger)
+
+
+class XASSpectrum(AbsorptionSpectrum):
+    """
+    X-ray Absorption Spectrum (XAS).
+    """
+
+    # ! implement `iri` and `rank` as part of `m_def = Section()`
+
+    xanes_spectrum = SubSection(
+        sub_section=AbsorptionSpectrum.m_def,
+        description="""
+        X-ray Absorption Near Edge Structure (XANES) spectrum.
+        """,
+        repeats=False,
+    )
+
+    exafs_spectrum = SubSection(
+        sub_section=AbsorptionSpectrum.m_def,
+        description="""
+        Extended X-ray Absorption Fine Structure (EXAFS) spectrum.
+        """,
+        repeats=False,
+    )
+
+    def __init__(
+        self, m_def: Section = None, m_context: Context = None, **kwargs
+    ) -> None:
+        super().__init__(m_def, m_context, **kwargs)
+        # Set the name of the section
+        self.name = self.m_def.name
+
+    def generate_from_contributions(self, logger: BoundLogger) -> None:
+        """
+        Generate the `value` of the XAS spectrum by concatenating the XANES and EXAFS contributions. It also concatenates
+        the `Energy` grid points of the XANES and EXAFS parts.
+
+        Args:
+            logger (BoundLogger): The logger to log messages.
+        """
+        # TODO check if this method is general enough
+        if self.xanes_spectrum is not None and self.exafs_spectrum is not None:
+            # Concatenate XANE and EXAFS `Energy` grid points
+            xanes_variables = get_variables(self.xanes_spectrum.variables, Energy)
+            exafs_variables = get_variables(self.exafs_spectrum.variables, Energy)
+            if len(xanes_variables) == 0 or len(exafs_variables) == 0:
+                logger.warning(
+                    'Could not extract the `Energy` grid points from XANES or EXAFS.'
+                )
+                return
+            xanes_energies = xanes_variables[0].points
+            exafs_energies = exafs_variables[0].points
+            if xanes_energies.max() > exafs_energies.min():
+                logger.warning(
+                    'The XANES `Energy` grid points are not below the EXAFS `Energy` grid points.'
+                )
+                return
+            self.variables = [
+                Energy(points=np.concatenate([xanes_energies, exafs_energies]))
+            ]
+            # Concatenate XANES and EXAFS `value` if they have the same shape ['n_energies']
+            try:
+                self.value = np.concatenate(
+                    [self.xanes_spectrum.value, self.exafs_spectrum.value]
+                )
+            except ValueError:
+                logger.warning(
+                    'The XANES and EXAFS `value` have different shapes. Could not concatenate the values.'
+                )
+
+    def normalize(self, archive, logger) -> None:
+        super().normalize(archive, logger)
+
+        if self.value is None:
+            logger.info(
+                'The `XASSpectrum.value` is not stored. We will attempt to obtain it by combining the XANES and EXAFS parts if these are present.'
+            )
+            self.generate_from_contributions(logger)
+