diff --git a/docs/digital_analog_qc/rydberg-hea.md b/docs/digital_analog_qc/rydberg-hea.md
deleted file mode 100644
index 1530c54c..00000000
--- a/docs/digital_analog_qc/rydberg-hea.md
+++ /dev/null
@@ -1,116 +0,0 @@
-Qadence simplifies the execution of digital-analog
-workloads on neutral atom quantum computers where the local
-addressability is restricted.
-
-In this regime, which we will refer to as *semi-local addressing*,
-the full Hamiltonian of the qubit system realized with neutral
-atoms comprises the following terms:
-
-$$
-\mathcal{H} = \mathcal{H}_{\textrm{global}} + \mathcal{H}_{\textrm{int}} + \mathcal{H}_{\textrm{local}}
-$$
-
-The first two terms are the standard components of a neutral atom Hamiltonians
-and read as follows:
-
-$$
-\mathcal{H}_{\textrm{global}} = \frac{\Omega}{2}\sum_{i}^N \left(
-    \textrm{cos}(\phi)\sigma^x_i - \textrm{sin}(\phi)\sigma^y_i \right) -
-    \delta \sum_{i}^N \hat{n}_i \\
-\mathcal{H}_{\textrm{int}} = \sum_{i<j} \frac{C_6}{|R_i - R_j|^6} \hat{n}_i \hat{n}_j
-$$
-
-where $\Omega$ is the Rabi frequency, $\phi$ the global phase, $\delta$ the
-detuning which can all be time-dependent (here omitted for simplicity). The operator
-$\hat{n}_i = \frac{1+\sigma_i^z}{2}$ is the occupation operator for the Rydberg state.
-$R_i$ represents instead the spatial coordinates of the i-th qubit.
-
-The local addressability term reads instead:
-
-$$
-\mathcal{H}_{\textrm{local}}(w^{drv}, w^{det}) = \frac{\tilde{\Omega}}{2}\sum_{i}^N
-    w_i^{drv}\left(\textrm{cos}(\phi)\sigma^x_i - \textrm{sin}(\phi)\sigma^y_i \right) -
-    \Delta \sum_{i}^N w_i^{det} \hat{n}_i
-$$
-
-In this Hamiltonian, the local addressing pattern in both Rabi frequency and detuning
-is determined by the weights $w^{drv} = \{w_i^{drv}\}$ and $w^{det} = \{w_i^{det}\}$ respectively. These
-weights are assigned before starting the simulation and they should have a unit sum.
-Their action is to effectively modulate the amplitude of the local drive/detuning pulses given by
-$\tilde{\Omega}$ and $\Delta$ which are here considered time-independent for simplicity.
-
-Qadence implements the Hamiltonian above in two different flavors of increasing complexity described
-below.
-
-## Circuit constructor
-
-The `rydberg_hea` constructor routine allows to
-build a circuit instance implementing a basic version of the Hamiltonian
-evolution described above where both $\Delta$ and $\tilde{\Omega}$ coefficients
-are considered constants. Furthemore, **no global drive and detuning** are explicitly added
-to the Hamiltonian. Therefore, the final Hamiltonian generator of the circuit reads as follows:
-
-$$
-\mathcal{H} = \mathcal{H}_{\textrm{local}}(w^{drv}, w^{det}) + \mathcal{H}_{\textrm{int}}
-$$
-
-This implementation does not perform any checks on the weights normalization, thus
-making it not realistic. This implies that global drive and detuning can be retrieved
-by appropriately choosing the weights.
-
-You can easily create a Rydberg hardware efficient ansatz implementing multiple layers
-of the evolution generated by the local addressing Hamiltonian:
-
-$$
-\mathcal{H}_{evo} = \sum_j \mathcal{H}_{\textrm{local}}(w_{j}^{drv}, w_{j}^{det})
-$$
-
-Notice that in real-device implementation, one layer only is usually possible.
-
-```python exec="on" source="material-block" result="json" session="rydb-hea"
-import qadence as qd
-from qadence import rydberg_hea, rydberg_hea_layer
-
-n_qubits = 4
-n_layers = 2
-register = qd.Register.line(n_qubits)
-
-# ansatz constructor
-# the evolution time is parametrized for each layer of the evolution
-ansatz = rydberg_hea(
-    register,
-    n_layers=n_layers,  # number of subsequent layers of Hamiltonian evolution
-    addressable_detuning=True,  # make the local detuning weights w_i^{det} as variational parameters
-    addressable_drive=True, # make the local drive weights w_i^{drv} as variational parameters
-    tunable_phase=True, # make the phase \phi as a variational parameter
-)
-
-# alternatively, a single ansatz layer can also be created for
-# better flexibility
-
-# these can be variational parameters
-tevo_drive = 1.0  # evolution time for the locally addressed drive term
-tevo_det = 1.0 # evolution time for the locally addressed detuning term
-tevo_int = 1.0  # evolution time for the interaction term
-
-# these can be list of variational parameters
-weights_drive = [0.0, 0.25, 0.5, 0.25]
-weights_det = [0.0, 0.0, 0.5, 0.5]
-
-ansatz_layer = rydberg_hea_layer(
-    register,
-    tevo_det,
-    tevo_drive,
-    tevo_int,
-    detunings=weights_det,
-    drives=weights_drive,
-)
-```
-
-This circuit constructor is meant to be used with fully differentiable backends such as
-PyQTorch and mainly for quick experimentation with neutral atom compatible ansatze.
-
-## Usage with digital-analog emulation
-
-A full integration with the [emulated digital-analog framework](./analog-basics.md) for realistic
-simulations is *coming soon*.
diff --git a/docs/digital_analog_qc/semi-local-addressing.md b/docs/digital_analog_qc/semi-local-addressing.md
index 5edbb860..04e58402 100644
--- a/docs/digital_analog_qc/semi-local-addressing.md
+++ b/docs/digital_analog_qc/semi-local-addressing.md
@@ -10,23 +10,100 @@ $$
 
 as is described in detail in the [analog interface basics](analog-basics.md) documentation.
 
-The driving Hamiltonian term in priciple can model any local single-qubit rotation by addressing each qubit individually. However, on the current generation of neutral-atom devices full local addressing is not yet achievable. Fortunatelly, using devices called spatial light modulators (SLMs) it is possible to implement semi-local addressing where the pattern of targeted qubits is fixed during the execution of the quantum circuit. More formally, the addressing pattern appears as an additional term in the neutral-atom Hamiltonian:
+The driving Hamiltonian term in priciple can model any local single-qubit rotation by addressing each qubit individually. However, some neutral-atom devices offer restricted local addressability using devices called spatial light modulators (SLMs).
+
+We refer to this regime as *semi-local addressability*. In this regime, the individual qubit addressing is restricted to a pattern of targeted
+qubits which is kept fixed during the execution of the quantum circuit. More formally, the addressing pattern appears as an additional
+term in the neutral-atom Hamiltonian:
 
 $$
-\mathcal{H} = \mathcal{H}_{\rm drive} + \mathcal{H}_{\rm int} + \mathcal{H}_{\rm pattern},
+\mathcal{H} = \mathcal{H}_{\rm drive} + \mathcal{H}_{\rm int} + \mathcal{H}_{\rm local}
 $$
 
 where $\mathcal{H}_{\rm pattern}$ is given by
 
 $$
-\mathcal{H}_{\rm pattern} = \sum_{i=0}^{n-1}\left(-\Delta w_i^{\rm det} N_i + \Gamma w_i^{\rm drive} X_i\right).
+\mathcal{H}_{\rm local} = \sum_{i=0}^{n-1}\left(-\Delta w_i^{\rm det} \hat{n}_i + \Gamma w_i^{\rm drive} \hat{\sigma}^x_i\right).
 $$
 
 Here $\Delta$ specifies the maximal **negative** detuning that each qubit in the register can be exposed to. The weight $w_i^{\rm det}\in [0, 1]$ determines the actual value of detuning that $i$-th qubit feels and this way the detuning pattern is emulated. Similarly, for the amplitude pattern $\Gamma$ determines the maximal additional **positive** drive that acts on qubits. In this case the corresponding weights $w_i^{\rm drive}$ can vary in the interval $[0, 1]$.
 
 Using the detuning and amplitude patterns described above one can modify the behavior of a selected set of qubits, thus achieving semi-local addressing.
 
-## Creating semi-local addressing patterns
+Qadence implements semi-local addressing in two different
+flavors of increasing complexity: either as a circuit
+constructor or directly as a pattern added to the general
+evolution Hamiltonian described by the circuit.
+
+## Using circuit constructors
+
+The `rydberg_hea` constructor routine allows to
+build a circuit instance implementing a basic version of the Hamiltonian
+evolution described above where both $\Delta$ and $\tilde{\Omega}$ coefficients
+are considered constants. Furthemore, **no global drive and detuning** are explicitly added
+to the Hamiltonian. Therefore, the final Hamiltonian generator of the circuit reads as follows:
+
+$$
+\mathcal{H} = \mathcal{H}_{\rm local}(w^{\rm drive}, w^{\rm det}) + \mathcal{H}_{\textrm{int}}
+$$
+
+This implementation does not perform any checks on the weights normalization, thus
+making it not realistic. This implies that global drive and detuning can be retrieved
+by appropriately choosing the weights.
+
+You can easily create a Rydberg hardware efficient ansatz implementing multiple layers
+of the evolution generated by the local addressing Hamiltonian:
+
+$$
+\mathcal{H}_{\rm evo} = \sum_j \mathcal{H}_{\textrm{local}}(w_{j}^{\rm drive}, w_{j}^{\rm det})
+$$
+
+Notice that in real-device implementation, one layer only is usually achievable.
+
+```python exec="on" source="material-block" result="json" session="addr"
+import qadence as qd
+from qadence import rydberg_hea, rydberg_hea_layer
+
+n_qubits = 4
+n_layers = 2
+register = qd.Register.line(n_qubits)
+
+# ansatz constructor
+# the evolution time is parametrized for each layer of the evolution
+ansatz = rydberg_hea(
+    register,
+    n_layers=n_layers,  # number of subsequent layers of Hamiltonian evolution
+    addressable_detuning=True,  # make the local detuning weights w_i^{det} as variational parameters
+    addressable_drive=True, # make the local drive weights w_i^{drv} as variational parameters
+    tunable_phase=True, # make the phase \phi as a variational parameter
+)
+
+# alternatively, a single ansatz layer can also be created for
+# better flexibility
+
+# these can be variational parameters
+tevo_drive = 1.0  # evolution time for the locally addressed drive term
+tevo_det = 1.0 # evolution time for the locally addressed detuning term
+tevo_int = 1.0  # evolution time for the interaction term
+
+# these can be list of variational parameters
+weights_drive = [0.0, 0.25, 0.5, 0.25]
+weights_det = [0.0, 0.0, 0.5, 0.5]
+
+ansatz_layer = rydberg_hea_layer(
+    register,
+    tevo_det,
+    tevo_drive,
+    tevo_int,
+    detunings=weights_det,
+    drives=weights_drive,
+)
+```
+
+This circuit constructor is meant to be used with fully differentiable backends such as
+`pyqtorch` and mainly for quick experimentation with neutral atom compatible ansatze.
+
+## Using addressing patterns
 
 In Qadence semi-local addressing patterns can be created by either specifying fixed values for the weights of the qubits being addressed or defining them as trainable parameters that can be optimized later in some training loop. Semi-local addressing patterns can be defined with the `AddressingPattern` dataclass.
 
@@ -34,7 +111,7 @@ In Qadence semi-local addressing patterns can be created by either specifying fi
 
 With fixed weights, detuning/amplitude addressing patterns can be defined in the following way:
 
-```python exec="on" source="material-block" session="emu"
+```python exec="on" source="material-block" session="addr"
 import torch
 from qadence.analog import AddressingPattern
 
@@ -58,7 +135,8 @@ If only detuning or amplitude pattern is needed - the corresponding weights for
 The created addressing pattern can now be passed as an argument to any Qadence device class, or to the
 `IdealDevice` or `RealisticDevice` to make use of the pre-defined options in those devices,
 
-```python exec="on" source="material-block" session="emu"
+```python exec="on" source="material-block" session="addr"
+import torch
 from qadence import (
     AnalogRX,
     AnalogRY,
@@ -102,7 +180,7 @@ print(f"Expectation value on PyQ: \n{expval_pyq.flatten().detach()}\n")  # markd
 
 The same configuration can also be seamlessly used to create a model with the Pulser backend.
 
-```python exec="on" source="material-block" session="emu"
+```python exec="on" source="material-block" session="addr"
 model_pulser = QuantumModel(
     circuit=circ,
     observable=obs,
@@ -130,7 +208,7 @@ it will be added to all the blocks.
 Since both the maximum detuning/amplitude value of the addressing pattern and the corresponding weights can be
 user specified, they can be variationally used in some QML setting. This can be achieved by defining pattern weights as trainable `Parameter` instances or strings specifying weight names.
 
-```python exec="on" source="material-block" session="emu"
+```python exec="on" source="material-block" session="addr"
 n_qubits = 3
 reg = Register.line(n_qubits, spacing=8.0)
 
diff --git a/mkdocs.yml b/mkdocs.yml
index 5249db6b..8128691d 100644
--- a/mkdocs.yml
+++ b/mkdocs.yml
@@ -23,11 +23,9 @@ nav:
     - digital_analog_qc/index.md
     - Basic operations on neutral-atoms: digital_analog_qc/analog-basics.md
     - Fitting a simple function: digital_analog_qc/analog-qcl.md
-    - Semi-local addressing pattern: digital_analog_qc/semi-local-addressing.md
-    - Hardware efficient ansatz with restricted addressability: digital_analog_qc/rydberg-hea.md
+    - Restricted local addressability: digital_analog_qc/semi-local-addressing.md
     - Pulse-level programming with Pulser: digital_analog_qc/pulser-basic.md
     - Solve a QUBO problem: digital_analog_qc/analog-qubo.md
-    - Pulse-level programming with Pulser: digital_analog_qc/pulser-basic.md
     - CNOT with interacting qubits: digital_analog_qc/daqc-cnot.md
 
   - Variational quantum algorithms:
diff --git a/qadence/constructors/__init__.py b/qadence/constructors/__init__.py
index 734a8a01..53bd6cf4 100644
--- a/qadence/constructors/__init__.py
+++ b/qadence/constructors/__init__.py
@@ -23,6 +23,7 @@
 )
 
 from .rydberg_hea import rydberg_hea, rydberg_hea_layer
+from .rydberg_feature_maps import rydberg_feature_map, analog_feature_map, rydberg_tower_feature_map
 
 from .qft import qft
 
@@ -45,4 +46,7 @@
     "daqc_transform",
     "rydberg_hea",
     "rydberg_hea_layer",
+    "rydberg_feature_map",
+    "analog_feature_map",
+    "rydberg_tower_feature_map",
 ]
diff --git a/qadence/constructors/feature_maps.py b/qadence/constructors/feature_maps.py
index b0754b37..84b94e24 100644
--- a/qadence/constructors/feature_maps.py
+++ b/qadence/constructors/feature_maps.py
@@ -6,7 +6,7 @@
 from math import isclose, pi
 from typing import Union
 
-from sympy import Function, acos
+from sympy import Basic, Function, acos
 
 from qadence.blocks import AbstractBlock, KronBlock, chain, kron, tag
 from qadence.logger import get_logger
@@ -36,65 +36,12 @@ def _set_range(fm_type: BasisSet | type[Function] | str) -> tuple[float, float]:
 }
 
 
-def feature_map(
-    n_qubits: int,
-    support: tuple[int, ...] | None = None,
+def fm_parameter(
+    fm_type: BasisSet | type[Function] | str,
     param: Parameter | str = "phi",
-    op: RotationTypes = RX,
-    fm_type: BasisSet | type[Function] | str = BasisSet.FOURIER,
-    reupload_scaling: ReuploadScaling | Callable | str = ReuploadScaling.CONSTANT,
     feature_range: tuple[float, float] | None = None,
     target_range: tuple[float, float] | None = None,
-    multiplier: Parameter | TParameter | None = None,
-) -> KronBlock:
-    """Construct a feature map of a given type.
-
-    Arguments:
-        n_qubits: Number of qubits the feature map covers. Results in `support=range(n_qubits)`.
-        support: Puts one feature-encoding rotation gate on every qubit in `support`. n_qubits in
-            this case specifies the total overall qubits of the circuit, which may be wider than the
-            support itself, but not narrower.
-        param: Parameter of the feature map; you can pass a string or Parameter;
-            it will be set as non-trainable (FeatureParameter) regardless.
-        op: Rotation operation of the feature map; choose from RX, RY, RZ or PHASE.
-        fm_type: Basis set for data encoding; choose from `BasisSet.FOURIER` for Fourier
-            encoding, or `BasisSet.CHEBYSHEV` for Chebyshev polynomials of the first kind.
-        reupload_scaling: how the feature map scales the data that is re-uploaded for each qubit.
-            choose from `ReuploadScaling` enumeration or provide your own function with a single
-            int as input and int or float as output.
-        feature_range: range of data that the input data is assumed to come from.
-        target_range: range of data the data encoder assumes as the natural range. For example,
-            in Chebyshev polynomials it is (-1, 1), while for Fourier it may be chosen as (0, 2*pi).
-        multiplier: overall multiplier; this is useful for reuploading the feature map serially with
-            different scalings; can be a number or parameter/expression.
-
-    Example:
-    ```python exec="on" source="material-block" result="json"
-    from qadence import feature_map, BasisSet, ReuploadScaling
-
-    fm = feature_map(3, fm_type=BasisSet.FOURIER)
-    print(f"{fm = }")
-
-    fm = feature_map(3, fm_type=BasisSet.CHEBYSHEV)
-    print(f"{fm = }")
-
-    fm = feature_map(3, fm_type=BasisSet.FOURIER, reupload_scaling = ReuploadScaling.TOWER)
-    print(f"{fm = }")
-    ```
-    """
-
-    # Process input
-    if support is None:
-        support = tuple(range(n_qubits))
-    elif len(support) != n_qubits:
-        raise ValueError("Wrong qubit support supplied")
-
-    if op not in ROTATIONS:
-        raise ValueError(
-            f"Operation {op} not supported. "
-            f"Please provide one from {[rot.__name__ for rot in ROTATIONS]}."
-        )
-
+) -> Parameter | Basic:
     # Backwards compatibility
     if fm_type in ("fourier", "chebyshev", "tower"):
         logger.warning(
@@ -108,7 +55,6 @@ def feature_map(
             fm_type = BasisSet.CHEBYSHEV
         elif fm_type == "tower":
             fm_type = BasisSet.CHEBYSHEV
-            reupload_scaling = ReuploadScaling.TOWER
 
     if isinstance(param, Parameter):
         fparam = param
@@ -144,8 +90,12 @@ def feature_map(
             "the given feature parameter with."
         )
 
-    basis_tag = fm_type.value if isinstance(fm_type, BasisSet) else str(fm_type)
+    return transformed_feature
+
 
+def fm_reupload_scaling_fn(
+    reupload_scaling: ReuploadScaling | Callable | str = ReuploadScaling.CONSTANT,
+) -> tuple[Callable, str]:
     # Set reupload scaling function
     if callable(reupload_scaling):
         rs_func = reupload_scaling
@@ -163,8 +113,82 @@ def feature_map(
         else:
             rs_tag = reupload_scaling
 
+    return rs_func, rs_tag
+
+
+def feature_map(
+    n_qubits: int,
+    support: tuple[int, ...] | None = None,
+    param: Parameter | str = "phi",
+    op: RotationTypes = RX,
+    fm_type: BasisSet | type[Function] | str = BasisSet.FOURIER,
+    reupload_scaling: ReuploadScaling | Callable | str = ReuploadScaling.CONSTANT,
+    feature_range: tuple[float, float] | None = None,
+    target_range: tuple[float, float] | None = None,
+    multiplier: Parameter | TParameter | None = None,
+) -> KronBlock:
+    """Construct a feature map of a given type.
+
+    Arguments:
+        n_qubits: Number of qubits the feature map covers. Results in `support=range(n_qubits)`.
+        support: Puts one feature-encoding rotation gate on every qubit in `support`. n_qubits in
+            this case specifies the total overall qubits of the circuit, which may be wider than the
+            support itself, but not narrower.
+        param: Parameter of the feature map; you can pass a string or Parameter;
+            it will be set as non-trainable (FeatureParameter) regardless.
+        op: Rotation operation of the feature map; choose from RX, RY, RZ or PHASE.
+        fm_type: Basis set for data encoding; choose from `BasisSet.FOURIER` for Fourier
+            encoding, or `BasisSet.CHEBYSHEV` for Chebyshev polynomials of the first kind.
+        reupload_scaling: how the feature map scales the data that is re-uploaded for each qubit.
+            choose from `ReuploadScaling` enumeration or provide your own function with a single
+            int as input and int or float as output.
+        feature_range: range of data that the input data is assumed to come from.
+        target_range: range of data the data encoder assumes as the natural range. For example,
+            in Chebyshev polynomials it is (-1, 1), while for Fourier it may be chosen as (0, 2*pi).
+        multiplier: overall multiplier; this is useful for reuploading the feature map serially with
+            different scalings; can be a number or parameter/expression.
+
+    Example:
+    ```python exec="on" source="material-block" result="json"
+    from qadence import feature_map, BasisSet, ReuploadScaling
+
+    fm = feature_map(3, fm_type=BasisSet.FOURIER)
+    print(f"{fm = }")
+
+    fm = feature_map(3, fm_type=BasisSet.CHEBYSHEV)
+    print(f"{fm = }")
+
+    fm = feature_map(3, fm_type=BasisSet.FOURIER, reupload_scaling = ReuploadScaling.TOWER)
+    print(f"{fm = }")
+    ```
+    """
+
+    # Process input
+    if support is None:
+        support = tuple(range(n_qubits))
+    elif len(support) != n_qubits:
+        raise ValueError("Wrong qubit support supplied")
+
+    if op not in ROTATIONS:
+        raise ValueError(
+            f"Operation {op} not supported. "
+            f"Please provide one from {[rot.__name__ for rot in ROTATIONS]}."
+        )
+
+    transformed_feature = fm_parameter(
+        fm_type, param, feature_range=feature_range, target_range=target_range
+    )
+
+    # Backwards compatibility
+    if fm_type == "tower":
+        logger.warning("Forcing reupload scaling strategy to TOWER")
+        reupload_scaling = ReuploadScaling.TOWER
+
+    basis_tag = fm_type.value if isinstance(fm_type, BasisSet) else str(fm_type)
+    rs_func, rs_tag = fm_reupload_scaling_fn(reupload_scaling)
+
     # Set overall multiplier
-    multiplier = 1 if multiplier is None else multiplier
+    multiplier = 1 if multiplier is None else Parameter(multiplier)
 
     # Build feature map
     op_list = []
diff --git a/qadence/constructors/rydberg_feature_maps.py b/qadence/constructors/rydberg_feature_maps.py
new file mode 100644
index 00000000..4cc8c09a
--- /dev/null
+++ b/qadence/constructors/rydberg_feature_maps.py
@@ -0,0 +1,113 @@
+from __future__ import annotations
+
+from typing import Callable
+
+import numpy as np
+from sympy import Basic, Function
+
+from qadence.blocks import AnalogBlock, KronBlock, kron
+from qadence.constructors.feature_maps import fm_parameter
+from qadence.logger import get_logger
+from qadence.operations import AnalogRot, AnalogRX, AnalogRY, AnalogRZ
+from qadence.parameters import FeatureParameter, Parameter, VariationalParameter
+from qadence.types import BasisSet, ReuploadScaling, TParameter
+
+logger = get_logger(__file__)
+
+AnalogRotationTypes = [AnalogRX, AnalogRY, AnalogRZ]
+
+
+def rydberg_feature_map(
+    n_qubits: int,
+    param: str = "phi",
+    max_abs_detuning: float = 2 * np.pi * 10,
+    weights: list[float] | None = None,
+) -> KronBlock:
+    """Feature map using semi-local addressing patterns.
+
+    If not weights are specified, variational parameters are created
+    for the pattern
+
+    Args:
+        n_qubits (int): number of qubits
+        param: the name of the feature parameter
+        max_abs_detuning: maximum value of absolute detuning for each qubit. Defaulted at 10 MHz.
+        weights: a list of wegiths to assign to each qubit parameter in the feature map
+
+    Returns:
+        The block representing the feature map
+    """
+
+    tower_coeffs: list[float | Parameter]
+    tower_coeffs = (
+        [VariationalParameter(f"w_{param}_{i}") for i in range(n_qubits)]
+        if weights is None
+        else weights
+    )
+    tower_detuning = max_abs_detuning / (sum(tower_coeffs[i] for i in range(n_qubits)))
+
+    param = FeatureParameter(param)
+    duration = 1000 * param / tower_detuning
+    return kron(
+        AnalogRot(
+            duration=duration,
+            delta=-tower_detuning * tower_coeffs[i],
+            phase=0.0,
+            qubit_support=(i,),
+        )
+        for i in range(n_qubits)
+    )
+
+
+def rydberg_tower_feature_map(
+    n_qubits: int, param: str = "phi", max_abs_detuning: float = 2 * np.pi * 10
+) -> KronBlock:
+    weights = list(np.arange(1, n_qubits + 1))
+    return rydberg_feature_map(
+        n_qubits, param=param, max_abs_detuning=max_abs_detuning, weights=weights
+    )
+
+
+def analog_feature_map(
+    param: str = "phi",
+    op: Callable[[Parameter | Basic], AnalogBlock] = AnalogRX,
+    fm_type: BasisSet | type[Function] | str = BasisSet.FOURIER,
+    reupload_scaling: ReuploadScaling | Callable | str = ReuploadScaling.CONSTANT,
+    feature_range: tuple[float, float] | None = None,
+    target_range: tuple[float, float] | None = None,
+    multiplier: Parameter | TParameter | None = None,
+) -> AnalogBlock:
+    """Generate a fully analog feature map.
+
+    Args:
+        param: Parameter of the feature map; you can pass a string or Parameter;
+            it will be set as non-trainable (FeatureParameter) regardless.
+        op: type of operation. Choose among AnalogRX, AnalogRY, AnalogRZ or a custom
+            callable function returning an AnalogBlock instance
+        fm_type: Basis set for data encoding; choose from `BasisSet.FOURIER` for Fourier
+            encoding, or `BasisSet.CHEBYSHEV` for Chebyshev polynomials of the first kind.
+        reupload_scaling: how the feature map scales the data that is re-uploaded. Given that
+            this feature map uses analog rotations, the reuploading works by simply
+            adding additional operations with different scaling factors in the parameter.
+            Choose from `ReuploadScaling` enumeration, currently only CONSTANT works,
+            or provide your own function with the first argument being the given
+            operation `op` and the second argument the feature parameter
+        feature_range: range of data that the input data is assumed to come from.
+        target_range: range of data the data encoder assumes as the natural range. For example,
+            in Chebyshev polynomials it is (-1, 1), while for Fourier it may be chosen as (0, 2*pi).
+        multiplier: overall multiplier; this is useful for reuploading the feature map serially with
+            different scalings; can be a number or parameter/expression.
+    """
+    transformed_feature = fm_parameter(
+        fm_type, param, feature_range=feature_range, target_range=target_range
+    )
+    multiplier = 1.0 if multiplier is None else Parameter(multiplier)
+
+    if callable(reupload_scaling):
+        return reupload_scaling(op, multiplier * transformed_feature)  # type: ignore[no-any-return]
+    elif reupload_scaling == ReuploadScaling.CONSTANT:
+        return op(multiplier * transformed_feature)
+    # TODO: implement tower scaling by reuploading multiple times
+    # using different analog rotations
+    else:
+        raise NotImplementedError(f"Reupload scaling {str(reupload_scaling)} not implemented!")
diff --git a/tests/constructors/test_rydberg_hea.py b/tests/constructors/test_rydberg_hea.py
index aa74d0b1..8dc0c903 100644
--- a/tests/constructors/test_rydberg_hea.py
+++ b/tests/constructors/test_rydberg_hea.py
@@ -1,9 +1,11 @@
 from __future__ import annotations
 
+import numpy as np
 import pytest
+import torch
 
 import qadence as qd
-from qadence import rydberg_hea
+from qadence import analog_feature_map, rydberg_feature_map, rydberg_hea, rydberg_tower_feature_map
 
 
 @pytest.mark.parametrize("detunings", [True, False])
@@ -62,3 +64,49 @@ def test_rydberg_hea_differentiation() -> None:
     for p in model.parameters():
         if p.requires_grad:
             assert p.grad is not None
+
+
+@pytest.mark.parametrize("basis", [qd.BasisSet.FOURIER, qd.BasisSet.CHEBYSHEV])
+def test_analog_feature_map(basis: qd.BasisSet) -> None:
+    pname = "x"
+    mname = "mult"
+    fm = analog_feature_map(
+        param=pname, op=qd.AnalogRY, fm_type=basis, multiplier=qd.VariationalParameter(mname)
+    )
+    assert isinstance(fm, qd.ConstantAnalogRotation)
+    assert fm.parameters.phase == -np.pi / 2
+    assert fm.parameters.delta == 0.0
+
+    params = list(fm.parameters.alpha.free_symbols)
+    assert len(params) == 2
+    assert pname in params and mname in params
+
+
+@pytest.mark.parametrize("weights", [None, [1.0, 2.0, 3.0, 4.0]])
+def test_rydberg_feature_map(weights: list[float] | None) -> None:
+    n_qubits = 4
+
+    fm = rydberg_feature_map(n_qubits, param="x", weights=weights)
+    assert len(fm) == n_qubits
+    assert all([isinstance(b, qd.ConstantAnalogRotation) for b in fm.blocks])
+
+    circuit = qd.QuantumCircuit(n_qubits, fm)
+    observable = qd.total_magnetization(n_qubits)
+    model = qd.QuantumModel(circuit, observable=observable)
+
+    values = {"x": torch.rand(1)}
+    expval = model.expectation(values)
+    expval.backward()
+    for p in model.parameters():
+        if p.requires_grad:
+            assert p.grad is not None
+
+
+def test_rydberg_tower_feature_map() -> None:
+    n_qubits = 4
+
+    fm1 = rydberg_tower_feature_map(n_qubits, param="x")
+    fm2 = rydberg_feature_map(n_qubits, param="x", weights=[1.0, 2.0, 3.0, 4.0])
+
+    for b1, b2 in zip(fm1.blocks, fm2.blocks):
+        assert b1.parameters.alpha == b2.parameters.alpha  # type:ignore [attr-defined]