diff --git a/src/qibocal/protocols/characterization/couplers/coupler_qubit_spectroscopy.py b/src/qibocal/protocols/characterization/couplers/coupler_qubit_spectroscopy.py
index 2ed4de04b..bfc608f17 100644
--- a/src/qibocal/protocols/characterization/couplers/coupler_qubit_spectroscopy.py
+++ b/src/qibocal/protocols/characterization/couplers/coupler_qubit_spectroscopy.py
@@ -48,13 +48,16 @@ def _acquisition(
     sequence = PulseSequence()
     ro_pulses = {}
     qd_pulses = {}
+    offset = {}
     couplers = []
     for i, pair in enumerate(targets):
         ordered_pair = order_pair(pair, platform)
         measured_qubit = params.measured_qubits[i]
 
         qubit = platform.qubits[measured_qubit].name
-        couplers.append(platform.pairs[tuple(sorted(ordered_pair))].coupler)
+        offset[qubit] = platform.pairs[tuple(sorted(ordered_pair))].coupler.sweetspot
+        coupler = platform.pairs[tuple(sorted(ordered_pair))].coupler.name
+        couplers.append(coupler)
 
         ro_pulses[qubit] = platform.create_qubit_readout_pulse(
             qubit, start=params.drive_duration
@@ -84,9 +87,11 @@ def _acquisition(
         (
             delta_bias_flux_range,
             sweepers,
+            sequences,
         ) = resonator_flux_dependence.create_flux_pulse_sweepers(
             params, platform, couplers, sequence
         )
+        sequence = sequences[0]
     else:
         delta_bias_flux_range = np.arange(
             -params.bias_width / 2, params.bias_width / 2, params.bias_step
@@ -103,6 +108,7 @@ def _acquisition(
     data = CouplerSpectroscopyData(
         resonator_type=platform.resonator_type,
         flux_pulses=params.flux_pulses,
+        offset=offset,
     )
 
     for bias_sweeper in sweepers:
diff --git a/src/qibocal/protocols/characterization/couplers/coupler_resonator_spectroscopy.py b/src/qibocal/protocols/characterization/couplers/coupler_resonator_spectroscopy.py
index 438fbb218..880824578 100644
--- a/src/qibocal/protocols/characterization/couplers/coupler_resonator_spectroscopy.py
+++ b/src/qibocal/protocols/characterization/couplers/coupler_resonator_spectroscopy.py
@@ -43,13 +43,15 @@ def _acquisition(
 
     sequence = PulseSequence()
     ro_pulses = {}
+    offset = {}
     couplers = []
     for i, pair in enumerate(targets):
         ordered_pair = order_pair(pair, platform)
         measured_qubit = params.measured_qubits[i]
 
         qubit = platform.qubits[measured_qubit].name
-        coupler = platform.pairs[tuple(sorted(ordered_pair))].coupler
+        offset[qubit] = platform.pairs[tuple(sorted(ordered_pair))].coupler.sweetspot
+        coupler = platform.pairs[tuple(sorted(ordered_pair))].coupler.name
         couplers.append(coupler)
         # TODO: May measure both qubits on the pair
         ro_pulses[qubit] = platform.create_qubit_readout_pulse(qubit, start=0)
@@ -72,12 +74,12 @@ def _acquisition(
 
     if params.flux_pulses:
         # TODO: Add delay
-        (
-            delta_bias_flux_range,
-            sweepers,
-        ) = resonator_flux_dependence.create_flux_pulse_sweepers(
-            params, platform, couplers, sequence
+        (delta_bias_flux_range, sweepers, sequences) = (
+            resonator_flux_dependence.create_flux_pulse_sweepers(
+                params, platform, couplers, sequence
+            )
         )
+        sequence = sequences[0]
     else:
         delta_bias_flux_range = np.arange(
             -params.bias_width / 2, params.bias_width / 2, params.bias_step
@@ -94,6 +96,7 @@ def _acquisition(
     data = CouplerSpectroscopyData(
         resonator_type=platform.resonator_type,
         flux_pulses=params.flux_pulses,
+        offset=offset,
     )
 
     for bias_sweeper in sweepers:
diff --git a/src/qibocal/protocols/characterization/couplers/utils.py b/src/qibocal/protocols/characterization/couplers/utils.py
index a66d205d3..b2518a347 100644
--- a/src/qibocal/protocols/characterization/couplers/utils.py
+++ b/src/qibocal/protocols/characterization/couplers/utils.py
@@ -53,6 +53,8 @@ class CouplerSpectroscopyData(Data):
     """Resonator type."""
     flux_pulses: bool
     """True if sweeping flux pulses, False if sweeping bias."""
+    offset: dict[QubitId, float] = field(default_factory=dict)
+    """Qubit bias offset."""
     data: dict[QubitId, npt.NDArray[CouplerSpecType]] = field(default_factory=dict)
     """Raw data acquired."""
 
diff --git a/src/qibocal/protocols/characterization/flux_dependence/qubit_crosstalk.py b/src/qibocal/protocols/characterization/flux_dependence/qubit_crosstalk.py
index 67b793a06..4a91d561e 100644
--- a/src/qibocal/protocols/characterization/flux_dependence/qubit_crosstalk.py
+++ b/src/qibocal/protocols/characterization/flux_dependence/qubit_crosstalk.py
@@ -8,10 +8,14 @@
 from qibolab.pulses import PulseSequence
 from qibolab.qubits import QubitId
 from qibolab.sweeper import Parameter, Sweeper, SweeperType
+from scipy.optimize import curve_fit
 
-from qibocal.auto.operation import Results, Routine
+from qibocal import update
+from qibocal.auto.operation import Routine
+from qibocal.config import log
 
 from ..qubit_spectroscopy_ef import DEFAULT_ANHARMONICITY
+from ..utils import HZ_TO_GHZ, table_dict, table_html
 from . import utils
 from .qubit_flux_dependence import (
     QubitFluxData,
@@ -19,6 +23,7 @@
     QubitFluxResults,
     QubitFluxType,
 )
+from .qubit_flux_dependence import _fit as diagonal_fit
 from .resonator_flux_dependence import create_flux_pulse_sweepers
 
 
@@ -32,13 +37,24 @@ class QubitCrosstalkParameters(QubitFluxParameters):
     If given flux will be swept on the given qubits in a sequential fashion (n qubits will result to n different executions).
     Multiple qubits may be measured in each execution as specified by the ``qubits`` option in the runcard.
     """
+    # TODO: add voltage parameters to bias qubits off sweetspot (absolute)
 
 
 @dataclass
 class QubitCrosstalkData(QubitFluxData):
     """Crosstalk acquisition outputs when ``flux_qubits`` are given."""
 
-    data: dict[QubitId, dict[QubitId, npt.NDArray[QubitFluxType]]] = field(
+    sweetspot: dict[QubitId, float] = field(default_factory=dict)
+    """Sweetspot for each qubit."""
+    d: dict[QubitId, float] = field(default_factory=dict)
+    """Asymmetry for each qubit."""
+    voltage: dict[QubitId, float] = field(default_factory=dict)
+    """Voltage applied to each flux line."""
+    matrix_element: dict[QubitId, float] = field(default_factory=dict)
+    """Diagonal crosstalk matrix element."""
+    qubit_frequency: dict[QubitId, float] = field(default_factory=dict)
+    """Qubit frequency for each qubit."""
+    data: dict[tuple[QubitId, QubitId], npt.NDArray[QubitFluxType]] = field(
         default_factory=dict
     )
     """Raw data acquired for (qubit, qubit_flux) pairs saved in nested dictionaries."""
@@ -53,13 +69,45 @@ def register_qubit(self, qubit, flux_qubit, freq, bias, signal, phase):
         else:
             self.data[qubit, flux_qubit] = ar
 
+    @property
+    def diagonal(self) -> Optional[QubitFluxData]:
+        instance = QubitFluxData(
+            resonator_type=self.resonator_type,
+            flux_pulses=self.flux_pulses,
+            qubit_frequency=self.qubit_frequency,
+        )
+        for qubit in self.qubits:
+            try:
+                instance.data[qubit] = self.data[qubit, qubit]
+            except KeyError:
+                log.info(
+                    f"Diagonal acquisition not found for qubit {qubit}. Runcard values will be used to perform the off-diagonal fit."
+                )
+
+        if len(instance.data) > 0:
+            return instance
+        return QubitFluxData(
+            resonator_type=self.resonator_type,
+            flux_pulses=self.flux_pulses,
+            qubit_frequency=self.qubit_frequency,
+        )
+
 
 @dataclass
-class QubitCrosstalkResult(Results):
+class QubitCrosstalkResults(QubitFluxResults):
     """
     Qubit Crosstalk outputs.
     """
 
+    crosstalk_matrix: dict[QubitId, dict[QubitId, float]]
+    """Crosstalk matrix element."""
+    fitted_parameters: dict[tuple[QubitId, QubitId], dict]
+    """Fitted parameters for each couple target-flux qubit."""
+
+    def __contains__(self, key: QubitId):
+        """Checking if qubit is in crosstalk_matrix attribute."""
+        return key in self.crosstalk_matrix
+
 
 def _acquisition(
     params: QubitCrosstalkParameters,
@@ -68,14 +116,26 @@ def _acquisition(
 ) -> QubitCrosstalkData:
     """Data acquisition for Crosstalk Experiment."""
 
-    # taking advantage of multiplexing, apply the same set of gates to all qubits in parallel
+    # TODO: pass voltage as parameter
     sequence = PulseSequence()
     ro_pulses = {}
     qd_pulses = {}
+    sweetspots = {}
+    d = {}
+    voltage = {}
+    matrix_element = {}
+    qubit_frequency = {}
     for qubit in targets:
+        qubit_frequency[qubit] = platform.qubits[qubit].drive_frequency
         qd_pulses[qubit] = platform.create_qubit_drive_pulse(
             qubit, start=0, duration=params.drive_duration
         )
+        try:
+            d[qubit] = platform.qubits[qubit].asymmetry
+            matrix_element[qubit] = platform.qubits[qubit].crosstalk_matrix[qubit]
+            sweetspots[qubit] = voltage[qubit] = platform.qubits[qubit].sweetspot
+        except KeyError:
+            log.warning(f"Missing flux parameters for qubit {qubit}.")
 
         if params.transition == "02":
             if platform.qubits[qubit].anharmonicity:
@@ -104,17 +164,18 @@ def _acquisition(
     )
     # TODO : abstract common lines with qubit flux dep routine
     if params.flux_qubits is None:
-        flux_qubits = list(platform.qubits.keys())
+        flux_qubits = list(platform.qubits)
     else:
         flux_qubits = params.flux_qubits
     if params.flux_pulses:
-        delta_bias_flux_range, sweepers = create_flux_pulse_sweepers(
-            params, platform, targets, sequence
+        delta_bias_flux_range, sweepers, sequences = create_flux_pulse_sweepers(
+            params, platform, flux_qubits, sequence, crosstalk=True
         )
     else:
         delta_bias_flux_range = np.arange(
             -params.bias_width / 2, params.bias_width / 2, params.bias_step
         )
+        sequences = [sequence] * len(flux_qubits)
         sweepers = [
             Sweeper(
                 Parameter.bias,
@@ -125,7 +186,13 @@ def _acquisition(
             for flux_qubit in flux_qubits
         ]
     data = QubitCrosstalkData(
-        resonator_type=platform.resonator_type, flux_pulses=params.flux_pulses
+        resonator_type=platform.resonator_type,
+        sweetspot=sweetspots,
+        voltage=voltage,
+        matrix_element=matrix_element,
+        d=d,
+        qubit_frequency=qubit_frequency,
+        flux_pulses=params.flux_pulses,
     )
 
     options = ExecutionParameters(
@@ -134,7 +201,7 @@ def _acquisition(
         acquisition_type=AcquisitionType.INTEGRATION,
         averaging_mode=AveragingMode.CYCLIC,
     )
-    for flux_qubit, bias_sweeper in zip(flux_qubits, sweepers):
+    for flux_qubit, bias_sweeper, sequence in zip(flux_qubits, sweepers, sequences):
         results = platform.sweep(sequence, options, bias_sweeper, freq_sweeper)
         # retrieve the results for every qubit
         for qubit in targets:
@@ -155,14 +222,137 @@ def _acquisition(
     return data
 
 
-def _fit(data: QubitCrosstalkData) -> QubitCrosstalkResult:
-    return QubitCrosstalkResult()
+def _fit(data: QubitCrosstalkData) -> QubitCrosstalkResults:
+
+    crosstalk_matrix = {qubit: {} for qubit in data.qubit_frequency}
+    fitted_parameters = {}
+
+    voltage = {}
+    asymmetry = {}
+    sweetspot = {}
+    matrix_element = {}
+    qubit_frequency = {}
+
+    diagonal = diagonal_fit(data.diagonal)
+
+    for qubit in data.qubits:
+        condition = qubit in diagonal
+        voltage[qubit] = diagonal.sweetspot[qubit] if condition else data.voltage[qubit]
+        asymmetry[qubit] = diagonal.asymmetry[qubit] if condition else data.d[qubit]
+        sweetspot[qubit] = (
+            diagonal.sweetspot[qubit] if condition else data.sweetspot[qubit]
+        )
+        matrix_element[qubit] = (
+            diagonal.matrix_element[qubit] if condition else data.matrix_element[qubit]
+        )
+        qubit_frequency[qubit] = (
+            diagonal.frequency[qubit] if condition else data.matrix_element[qubit]
+        )
+
+    for target_flux_qubit, qubit_data in data.data.items():
+
+        if data.resonator_type == "3D":
+            frequencies, biases = utils.extract_min_feature(
+                qubit_data.freq,
+                qubit_data.bias,
+                qubit_data.signal,
+            )
+        else:
+            frequencies, biases = utils.extract_max_feature(
+                qubit_data.freq,
+                qubit_data.bias,
+                qubit_data.signal,
+            )
+        target_qubit, flux_qubit = target_flux_qubit
+
+        if target_qubit != flux_qubit:
+            # fit function needs to be defined here to pass correct parameters
+            # at runtime
+            def fit_function(x, crosstalk_element):
+                return utils.transmon_frequency(
+                    xi=voltage[target_qubit],
+                    xj=x,
+                    w_max=qubit_frequency[target_qubit] * HZ_TO_GHZ,
+                    d=asymmetry[target_qubit],
+                    sweetspot=sweetspot[target_qubit],
+                    matrix_element=matrix_element[target_qubit],
+                    crosstalk_element=crosstalk_element,
+                )
+
+            try:
+                popt, _ = curve_fit(
+                    fit_function,
+                    biases,
+                    frequencies * HZ_TO_GHZ,
+                    bounds=(-np.inf, np.inf),
+                )
+
+                fitted_parameters[target_qubit, flux_qubit] = dict(
+                    xi=voltage[target_qubit],
+                    w_max=qubit_frequency[target_qubit],
+                    d=asymmetry[target_qubit],
+                    sweetspot=sweetspot[target_qubit],
+                    matrix_element=matrix_element[target_qubit],
+                    crosstalk_element=float(popt),
+                )
+                crosstalk_matrix[target_qubit][flux_qubit] = float(popt)
+            except ValueError as e:
+                log.error(
+                    f"Off-diagonal flux fit failed for qubit {flux_qubit} due to {e}."
+                )
+
+        else:
+            fitted_parameters[target_qubit, flux_qubit] = diagonal.fitted_parameters[
+                target_qubit
+            ]
+            crosstalk_matrix[target_qubit][flux_qubit] = matrix_element[target_qubit]
+
+    return QubitCrosstalkResults(
+        frequency=qubit_frequency,
+        sweetspot=sweetspot,
+        asymmetry=asymmetry,
+        matrix_element=matrix_element,
+        crosstalk_matrix=crosstalk_matrix,
+        fitted_parameters=fitted_parameters,
+    )
 
 
-def _plot(data: QubitFluxData, fit: QubitFluxResults, target: QubitId):
+def _plot(data: QubitCrosstalkData, fit: QubitCrosstalkResults, target: QubitId):
     """Plotting function for Crosstalk Experiment."""
-    return utils.flux_crosstalk_plot(data, target)
+    figures, fitting_report = utils.flux_crosstalk_plot(
+        data, target, fit, fit_function=utils.transmon_frequency
+    )
+    if fit is not None:
+        labels = ["Sweetspot [V]", "Qubit Frequency at Sweetspot [Hz]", "Asymmetry d"]
+        values = [
+            np.round(fit.sweetspot[target], 4),
+            np.round(fit.frequency[target], 4),
+            np.round(fit.asymmetry[target], 4),
+        ]
+        for flux_qubit in fit.crosstalk_matrix[target]:
+            if flux_qubit != target:
+                labels.append(f"Crosstalk with qubit {flux_qubit}")
+            else:
+                labels.append(f"Flux dependence")
+            values.append(np.round(fit.crosstalk_matrix[target][flux_qubit], 4))
+        fitting_report = table_html(
+            table_dict(
+                target,
+                labels,
+                values,
+            )
+        )
+    return figures, fitting_report
+
+
+def _update(results: QubitCrosstalkResults, platform: Platform, qubit: QubitId):
+    """Update crosstalk matrix."""
+    update.drive_frequency(results.frequency[qubit], platform, qubit)
+    update.sweetspot(results.sweetspot[qubit], platform, qubit)
+    update.asymmetry(results.asymmetry[qubit], platform, qubit)
+    for flux_qubit, element in results.crosstalk_matrix[qubit].items():
+        update.crosstalk_matrix(element, platform, qubit, flux_qubit)
 
 
-qubit_crosstalk = Routine(_acquisition, _fit, _plot)
+qubit_crosstalk = Routine(_acquisition, _fit, _plot, _update)
 """Qubit crosstalk Routine object"""
diff --git a/src/qibocal/protocols/characterization/flux_dependence/qubit_flux_dependence.py b/src/qibocal/protocols/characterization/flux_dependence/qubit_flux_dependence.py
index 9f523b73a..0f53a0cf0 100644
--- a/src/qibocal/protocols/characterization/flux_dependence/qubit_flux_dependence.py
+++ b/src/qibocal/protocols/characterization/flux_dependence/qubit_flux_dependence.py
@@ -75,6 +75,9 @@ class QubitFluxData(Data):
     qubit_frequency: dict[QubitId, float] = field(default_factory=dict)
     """Qubit frequencies."""
 
+    offset: dict[QubitId, float] = field(default_factory=dict)
+    """Qubit bias offset."""
+
     data: dict[QubitId, npt.NDArray[QubitFluxType]] = field(default_factory=dict)
     """Raw data acquired."""
 
@@ -97,11 +100,13 @@ def _acquisition(
     ro_pulses = {}
     qd_pulses = {}
     qubit_frequency = {}
+    offset = {}
     for qubit in targets:
         qd_pulses[qubit] = platform.create_qubit_drive_pulse(
             qubit, start=0, duration=params.drive_duration
         )
         qubit_frequency[qubit] = platform.qubits[qubit].drive_frequency
+        offset[qubit] = platform.qubits[qubit].sweetspot
 
         if params.transition == "02":
             if platform.qubits[qubit].anharmonicity:
@@ -129,12 +134,12 @@ def _acquisition(
         type=SweeperType.OFFSET,
     )
     if params.flux_pulses:
-        (
-            delta_bias_flux_range,
-            sweepers,
-        ) = resonator_flux_dependence.create_flux_pulse_sweepers(
-            params, platform, targets, sequence
+        (delta_bias_flux_range, sweepers, sequences) = (
+            resonator_flux_dependence.create_flux_pulse_sweepers(
+                params, platform, targets, sequence
+            )
         )
+        sequence = sequences[0]
     else:
         delta_bias_flux_range = np.arange(
             -params.bias_width / 2, params.bias_width / 2, params.bias_step
@@ -151,8 +156,8 @@ def _acquisition(
         resonator_type=platform.resonator_type,
         flux_pulses=params.flux_pulses,
         qubit_frequency=qubit_frequency,
+        offset=offset,
     )
-
     options = ExecutionParameters(
         nshots=params.nshots,
         relaxation_time=params.relaxation_time,
@@ -200,7 +205,7 @@ def _fit(data: QubitFluxData) -> QubitFluxResults:
 
         try:
             popt = curve_fit(
-                utils.transmon_frequency,
+                utils.transmon_frequency_diagonal,
                 biases,
                 frequencies * HZ_TO_GHZ,
                 bounds=utils.qubit_flux_dependence_fit_bounds(
@@ -233,7 +238,10 @@ def _fit(data: QubitFluxData) -> QubitFluxResults:
 def _plot(data: QubitFluxData, fit: QubitFluxResults, target: QubitId):
     """Plotting function for QubitFlux Experiment."""
     figures = utils.flux_dependence_plot(
-        data, fit, target, fit_function=utils.transmon_frequency
+        data,
+        fit,
+        target,
+        fit_function=utils.transmon_frequency_diagonal,
     )
     if data.flux_pulses:
         bias_flux_unit = "a.u."
@@ -247,7 +255,7 @@ def _plot(data: QubitFluxData, fit: QubitFluxResults, target: QubitId):
                     f"Sweetspot [{bias_flux_unit}]",
                     "Qubit Frequency at Sweetspot [Hz]",
                     "Asymmetry d",
-                    "V_ii [V]",
+                    "Flux dependence",
                 ],
                 [
                     np.round(fit.sweetspot[target], 4),
@@ -265,6 +273,7 @@ def _update(results: QubitFluxResults, platform: Platform, qubit: QubitId):
     update.drive_frequency(results.frequency[qubit], platform, qubit)
     update.sweetspot(results.sweetspot[qubit], platform, qubit)
     update.asymmetry(results.asymmetry[qubit], platform, qubit)
+    update.crosstalk_matrix(results.matrix_element[qubit], platform, qubit, qubit)
 
 
 qubit_flux = Routine(_acquisition, _fit, _plot, _update)
diff --git a/src/qibocal/protocols/characterization/flux_dependence/qubit_flux_tracking.py b/src/qibocal/protocols/characterization/flux_dependence/qubit_flux_tracking.py
index 0b0fc5d43..dc0f9160a 100644
--- a/src/qibocal/protocols/characterization/flux_dependence/qubit_flux_tracking.py
+++ b/src/qibocal/protocols/characterization/flux_dependence/qubit_flux_tracking.py
@@ -102,7 +102,7 @@ def _acquisition(
     for bias in delta_bias_range:
         for qubit in targets:
             try:
-                freq_resonator = utils.transmon_readout_frequency(
+                freq_resonator = utils.transmon_readout_frequency_diagonal(
                     bias,
                     platform.qubits[qubit].drive_frequency,
                     platform.qubits[qubit].asymmetry,
diff --git a/src/qibocal/protocols/characterization/flux_dependence/resonator_crosstalk.py b/src/qibocal/protocols/characterization/flux_dependence/resonator_crosstalk.py
index eb8c480cb..9bb697a85 100644
--- a/src/qibocal/protocols/characterization/flux_dependence/resonator_crosstalk.py
+++ b/src/qibocal/protocols/characterization/flux_dependence/resonator_crosstalk.py
@@ -8,16 +8,22 @@
 from qibolab.pulses import PulseSequence
 from qibolab.qubits import QubitId
 from qibolab.sweeper import Parameter, Sweeper, SweeperType
+from scipy.optimize import curve_fit
 
-from qibocal.auto.operation import Results, Routine
-from qibocal.protocols.characterization.flux_dependence.resonator_flux_dependence import (
+from qibocal import update
+from qibocal.auto.operation import Routine
+from qibocal.config import log
+
+from ..utils import HZ_TO_GHZ, table_dict, table_html
+from . import utils
+from .resonator_flux_dependence import (
     ResFluxType,
     ResonatorFluxData,
     ResonatorFluxParameters,
-    create_flux_pulse_sweepers,
+    ResonatorFluxResults,
 )
-
-from . import utils
+from .resonator_flux_dependence import _fit as diagonal_fit
+from .resonator_flux_dependence import create_flux_pulse_sweepers
 
 
 @dataclass
@@ -33,14 +39,35 @@ class ResCrosstalkParameters(ResonatorFluxParameters):
 
 
 @dataclass
-class ResCrosstalkResults(Results):
-    """Empty fitting outputs for cross talk because fitting is not implemented in this case."""
+class ResCrosstalkResults(ResonatorFluxResults):
+    """ResCrosstalk outputs."""
+
+    crosstalk_matrix: dict[QubitId, dict[QubitId, float]] = field(default_factory=dict)
+    """Crosstalk matrix element."""
+    fitted_parameters: dict[tuple[QubitId, QubitId], dict] = field(default_factory=dict)
+    """Fitted parameters for each couple target-flux qubit."""
+
+    def __contains__(self, key: QubitId):
+        """Checking if qubit is in crosstalk_matrix attribute."""
+        return key in self.crosstalk_matrix
 
 
 @dataclass
 class ResCrosstalkData(ResonatorFluxData):
-    """QubitFlux acquisition outputs when ``flux_qubits`` are given."""
+    """ResFlux acquisition outputs when ``flux_qubits`` are given."""
 
+    sweetspot: dict[QubitId, float] = field(default_factory=dict)
+    """Sweetspot for each qubit."""
+    asymmetry: dict[QubitId, float] = field(default_factory=dict)
+    """Asymmetry for each qubit."""
+    coupling: dict[QubitId, float] = field(default_factory=dict)
+    """Coupling parameter g for each qubit."""
+    voltage: dict[QubitId, float] = field(default_factory=dict)
+    """Voltage provided to each qubit."""
+    resonator_frequency: dict[QubitId, float] = field(default_factory=dict)
+    """Readout resonator frequency for each qubit."""
+    matrix_element: dict[QubitId, float] = field(default_factory=dict)
+    """Diagonal crosstalk matrix element."""
     data: dict[tuple[QubitId, QubitId], npt.NDArray[ResFluxType]] = field(
         default_factory=dict
     )
@@ -56,10 +83,34 @@ def register_qubit(self, qubit, flux_qubit, freq, bias, signal, phase):
         else:
             self.data[qubit, flux_qubit] = ar
 
+    @property
+    def diagonal(self) -> Optional[ResonatorFluxData]:
+        instance = ResonatorFluxData(
+            resonator_type=self.resonator_type,
+            flux_pulses=self.flux_pulses,
+            qubit_frequency=self.qubit_frequency,
+            bare_resonator_frequency=self.bare_resonator_frequency,
+        )
+        for qubit in self.qubits:
+            try:
+                instance.data[qubit] = self.data[qubit, qubit]
+            except KeyError:
+                log.info(
+                    f"Diagonal acquisition not found for qubit {qubit}. Runcard values will be used to perform the off-diagonal fit."
+                )
+
+        if len(instance.data) > 0:
+            return instance
+        return ResonatorFluxData(
+            resonator_type=self.resonator_type,
+            flux_pulses=self.flux_pulses,
+            qubit_frequency=self.qubit_frequency,
+        )
+
 
 def _acquisition(
     params: ResCrosstalkParameters, platform: Platform, targets: list[QubitId]
-) -> ResonatorFluxData:
+) -> ResCrosstalkData:
     """Data acquisition for ResonatorFlux experiment."""
     # create a sequence of pulses for the experiment:
     # MZ
@@ -68,12 +119,27 @@ def _acquisition(
     sequence = PulseSequence()
     ro_pulses = {}
     bare_resonator_frequency = {}
+    resonator_frequency = {}
     qubit_frequency = {}
+    sweetspots = {}
+    asymmetry = {}
+    coupling = {}
+    voltage = {}
+    matrix_element = {}
     for qubit in targets:
+        try:
+            sweetspots[qubit] = voltage[qubit] = platform.qubits[qubit].sweetspot
+            asymmetry[qubit] = platform.qubits[qubit].asymmetry
+            coupling[qubit] = platform.qubits[qubit].g
+            matrix_element[qubit] = platform.qubits[qubit].crosstalk_matrix[qubit]
+        except KeyError:
+            log.warning(f"Missing flux parameters for qubit {qubit}.")
+
         bare_resonator_frequency[qubit] = platform.qubits[
             qubit
         ].bare_resonator_frequency
         qubit_frequency[qubit] = platform.qubits[qubit].drive_frequency
+        resonator_frequency[qubit] = platform.qubits[qubit].readout_frequency
 
         ro_pulses[qubit] = platform.create_qubit_readout_pulse(qubit, start=0)
         sequence.add(ro_pulses[qubit])
@@ -90,18 +156,19 @@ def _acquisition(
     )
 
     if params.flux_qubits is None:
-        flux_qubits = list(platform.qubits.keys())
+        flux_qubits = list(platform.qubits)
 
     else:
         flux_qubits = params.flux_qubits
     if params.flux_pulses:
-        delta_bias_flux_range, sweepers = create_flux_pulse_sweepers(
-            params, platform, targets, sequence
+        delta_bias_flux_range, sweepers, sequences = create_flux_pulse_sweepers(
+            params, platform, flux_qubits, sequence, crosstalk=True
         )
     else:
         delta_bias_flux_range = np.arange(
             -params.bias_width / 2, params.bias_width / 2, params.bias_step
         )
+        sequences = [sequence] * len(flux_qubits)
         sweepers = [
             Sweeper(
                 Parameter.bias,
@@ -116,6 +183,12 @@ def _acquisition(
         resonator_type=platform.resonator_type,
         flux_pulses=params.flux_pulses,
         qubit_frequency=qubit_frequency,
+        resonator_frequency=resonator_frequency,
+        sweetspot=sweetspots,
+        voltage=voltage,
+        matrix_element=matrix_element,
+        asymmetry=asymmetry,
+        coupling=coupling,
         bare_resonator_frequency=bare_resonator_frequency,
     )
     options = ExecutionParameters(
@@ -124,7 +197,7 @@ def _acquisition(
         acquisition_type=AcquisitionType.INTEGRATION,
         averaging_mode=AveragingMode.CYCLIC,
     )
-    for flux_qubit, bias_sweeper in zip(flux_qubits, sweepers):
+    for flux_qubit, bias_sweeper, sequence in zip(flux_qubits, sweepers, sequences):
         results = platform.sweep(sequence, options, bias_sweeper, freq_sweeper)
         # retrieve the results for every qubit
         for qubit in targets:
@@ -141,18 +214,165 @@ def _acquisition(
                 freq=delta_frequency_range + ro_pulses[qubit].frequency,
                 bias=delta_bias_flux_range + sweetspot,
             )
-
     return data
 
 
 def _fit(data: ResCrosstalkData) -> ResCrosstalkResults:
-    return ResCrosstalkResults()
+    crosstalk_matrix = {qubit: {} for qubit in data.qubit_frequency}
+    fitted_parameters = {}
+    diagonal = diagonal_fit(data.diagonal)
+
+    voltage = {}
+    sweetspot = {}
+    asymmetry = {}
+    coupling = {}
+    matrix_element = {}
+    qubit_frequency = {}
+    bare_resonator_frequency = {}
+    resonator_frequency = {}
+
+    for qubit in data.qubits:
+        condition = qubit in diagonal
+        voltage[qubit] = diagonal.sweetspot[qubit] if condition else data.voltage[qubit]
+        sweetspot[qubit] = (
+            diagonal.sweetspot[qubit] if condition else data.sweetspot[qubit]
+        )
+        asymmetry[qubit] = (
+            diagonal.asymmetry[qubit] if condition else data.asymmetry[qubit]
+        )
+        coupling[qubit] = (
+            diagonal.coupling[qubit] if condition else data.coupling[qubit]
+        )
+        matrix_element[qubit] = (
+            diagonal.matrix_element[qubit] if condition else data.matrix_element[qubit]
+        )
+        qubit_frequency[qubit] = (
+            diagonal.drive_frequency[qubit]
+            if condition
+            else data.qubit_frequency[qubit]
+        )
+        bare_resonator_frequency[qubit] = (
+            diagonal.bare_frequency[qubit]
+            if condition
+            else data.bare_resonator_frequency[qubit]
+        )
+        resonator_frequency[qubit] = (
+            diagonal.frequency[qubit] if condition else data.resonator_frequency[qubit]
+        )
+
+    for target_flux_qubit, qubit_data in data.data.items():
+        target_qubit, flux_qubit = target_flux_qubit
+
+        if data.resonator_type == "3D":
+            frequencies, biases = utils.extract_max_feature(
+                qubit_data.freq,
+                qubit_data.bias,
+                qubit_data.signal,
+            )
+        else:
+            frequencies, biases = utils.extract_min_feature(
+                qubit_data.freq,
+                qubit_data.bias,
+                qubit_data.signal,
+            )
+
+        if target_qubit != flux_qubit:
+            # fit function needs to be defined here to pass correct parameters
+            # at runtime
+            def fit_function(x, crosstalk_element):
+                return utils.transmon_readout_frequency(
+                    xi=voltage[target_qubit],
+                    xj=x,
+                    w_max=qubit_frequency[target_qubit],
+                    d=asymmetry[target_qubit],
+                    sweetspot=sweetspot[target_qubit],
+                    matrix_element=matrix_element[target_qubit],
+                    g=coupling[target_qubit],
+                    resonator_freq=bare_resonator_frequency[target_qubit],
+                    crosstalk_element=crosstalk_element,
+                )
+
+            try:
+                popt, _ = curve_fit(
+                    fit_function, biases, frequencies * HZ_TO_GHZ, bounds=(-1e-1, 1e-1)
+                )
+                fitted_parameters[target_qubit, flux_qubit] = dict(
+                    xi=voltage[target_qubit],
+                    w_max=qubit_frequency[target_qubit],
+                    d=asymmetry[target_qubit],
+                    sweetspot=sweetspot[target_qubit],
+                    matrix_element=matrix_element[target_qubit],
+                    g=coupling[target_qubit],
+                    resonator_freq=bare_resonator_frequency[target_qubit],
+                    crosstalk_element=float(popt),
+                )
+                crosstalk_matrix[target_qubit][flux_qubit] = float(popt)
+            except ValueError as e:
+                log.error(
+                    f"Off-diagonal flux fit failed for qubit {flux_qubit} due to {e}."
+                )
+        else:
+            fitted_parameters[target_qubit, flux_qubit] = diagonal.fitted_parameters[
+                target_qubit
+            ]
+            crosstalk_matrix[target_qubit][flux_qubit] = matrix_element[target_qubit]
+
+    return ResCrosstalkResults(
+        frequency=resonator_frequency,
+        sweetspot=sweetspot,
+        asymmetry=asymmetry,
+        bare_frequency=bare_resonator_frequency,
+        drive_frequency=qubit_frequency,
+        coupling=coupling,
+        crosstalk_matrix=crosstalk_matrix,
+        fitted_parameters=fitted_parameters,
+    )
 
 
 def _plot(data: ResCrosstalkData, fit: ResCrosstalkResults, target: QubitId):
     """Plotting function for ResonatorFlux Experiment."""
-    return utils.flux_crosstalk_plot(data, target)
+    figures, fitting_report = utils.flux_crosstalk_plot(
+        data, target, fit, fit_function=utils.transmon_readout_frequency
+    )
+    if fit is not None:
+        labels = [
+            "Sweetspot [V]",
+            "Resonator Frequency at Sweetspot [Hz]",
+            "Asymmetry d",
+            "Coupling g",
+            "Bare Resonator Frequency [Hz]",
+            "Qubit Frequency [Hz]",
+        ]
+        values = [
+            np.round(fit.sweetspot[target], 4),
+            np.round(fit.frequency[target], 4),
+            np.round(fit.asymmetry[target], 4),
+            np.round(fit.coupling[target], 4),
+            np.round(fit.bare_frequency[target], 4),
+            np.round(fit.drive_frequency[target], 4),
+        ]
+        for flux_qubit in fit.crosstalk_matrix[target]:
+            if flux_qubit != target:
+                labels.append(f"Crosstalk with qubit {flux_qubit}")
+            else:
+                labels.append(f"Flux dependence")
+            values.append(np.round(fit.crosstalk_matrix[target][flux_qubit], 4))
+
+        fitting_report = table_html(
+            table_dict(
+                target,
+                labels,
+                values,
+            )
+        )
+    return figures, fitting_report
+
+
+def _update(results: ResCrosstalkResults, platform: Platform, qubit: QubitId):
+    """Update crosstalk matrix."""
+    for flux_qubit, element in results.crosstalk_matrix[qubit].items():
+        update.crosstalk_matrix(element, platform, qubit, flux_qubit)
 
 
-resonator_crosstalk = Routine(_acquisition, _fit, _plot)
+resonator_crosstalk = Routine(_acquisition, _fit, _plot, _update)
 """Resonator crosstalk Routine object"""
diff --git a/src/qibocal/protocols/characterization/flux_dependence/resonator_flux_dependence.py b/src/qibocal/protocols/characterization/flux_dependence/resonator_flux_dependence.py
index 1de96656a..6cdf7fb14 100644
--- a/src/qibocal/protocols/characterization/flux_dependence/resonator_flux_dependence.py
+++ b/src/qibocal/protocols/characterization/flux_dependence/resonator_flux_dependence.py
@@ -1,10 +1,10 @@
+from copy import deepcopy
 from dataclasses import dataclass, field
 from typing import List, Optional, Union
 
 import numpy as np
 import numpy.typing as npt
 from qibolab import AcquisitionType, AveragingMode, ExecutionParameters
-from qibolab.couplers import Coupler
 from qibolab.platform import Platform
 from qibolab.pulses import PulseSequence
 from qibolab.qubits import QubitId
@@ -12,7 +12,7 @@
 from scipy.optimize import curve_fit
 
 from qibocal import update
-from qibocal.auto.operation import Data, Parameters, Qubits, Results, Routine
+from qibocal.auto.operation import Data, Parameters, Results, Routine
 from qibocal.config import log
 
 from ..utils import GHZ_TO_HZ, HZ_TO_GHZ, table_dict, table_html
@@ -99,21 +99,21 @@ def flux_pulses(self):
 class ResonatorFluxResults(Results):
     """ResonatoFlux outputs."""
 
-    frequency: dict[QubitId, float]
+    frequency: dict[QubitId, float] = field(default_factory=dict)
     """Readout frequency for each qubit."""
-    sweetspot: dict[QubitId, float]
+    sweetspot: dict[QubitId, float] = field(default_factory=dict)
     """Sweetspot for each qubit."""
-    asymmetry: dict[QubitId, float]
+    asymmetry: dict[QubitId, float] = field(default_factory=dict)
     """Asymmetry between junctions."""
-    bare_frequency: dict[QubitId, float]
+    bare_frequency: dict[QubitId, float] = field(default_factory=dict)
     """Resonator bare frequency."""
-    drive_frequency: dict[QubitId, float]
+    drive_frequency: dict[QubitId, float] = field(default_factory=dict)
     """Qubit frequency at sweetspot."""
-    fitted_parameters: dict[QubitId, dict[str, float]]
+    fitted_parameters: dict[QubitId, dict[str, float]] = field(default_factory=dict)
     """Raw fitting output."""
-    coupling: dict[QubitId, float]
+    coupling: dict[QubitId, float] = field(default_factory=dict)
     """Qubit-resonator coupling."""
-    matrix_element: dict[QubitId, float]
+    matrix_element: dict[QubitId, float] = field(default_factory=dict)
     """C_ii coefficient."""
 
 
@@ -134,13 +134,12 @@ class ResonatorFluxData(Data):
 
     resonator_type: str
     """Resonator type."""
-
     flux_pulses: bool
     """True if sweeping flux pulses, False if sweeping bias."""
-
     qubit_frequency: dict[QubitId, float] = field(default_factory=dict)
     """Qubit frequencies."""
-
+    offset: dict[QubitId, float] = field(default_factory=dict)
+    """Qubit bias offset."""
     bare_resonator_frequency: dict[QubitId, int] = field(default_factory=dict)
     """Qubit bare resonator frequency power provided by the user."""
 
@@ -157,9 +156,10 @@ def register_qubit(self, qubit, freq, bias, signal, phase):
 def create_flux_pulse_sweepers(
     params: ResonatorFluxParameters,
     platform: Platform,
-    qubits: Qubits,
+    qubits: list[QubitId],
     sequence: PulseSequence,
-) -> tuple[np.ndarray, list[Sweeper]]:
+    crosstalk: bool = False,
+) -> tuple[np.ndarray, list[Sweeper], list[PulseSequence]]:
     """Create a list of sweepers containing flux pulses.
 
     Args:
@@ -167,8 +167,10 @@ def create_flux_pulse_sweepers(
         platform (Platform): platform on which to run the experiment.
         qubits (Qubits): qubits on which to run the experiment.
         sequence (PulseSequence): pulse sequence of the experiment (updated with flux pulses).
+        crosstalk (bool): if True it will split amplitude sweepers (necessary for crosstalk protocol)
     """
     qf_pulses = {}
+    sequences = [deepcopy(sequence) for _ in range(len(qubits))]
     for i, qubit in enumerate(qubits):
         if isinstance(params.flux_amplitude_start, list):
             flux_amplitude_start = params.flux_amplitude_start[i]
@@ -183,37 +185,46 @@ def create_flux_pulse_sweepers(
             flux_amplitude_end,
             flux_amplitude_step,
         )
-        if isinstance(qubit, Coupler):
+
+        if qubit not in platform.qubits:
             # FIXME: Missmatch with create_coupler_pulse and create_qubit_flux_pulse
             pulse = platform.create_coupler_pulse(
-                qubit.name,
+                qubit,
                 start=0,
                 duration=sequence.duration,
                 amplitude=1,
             )
-            qubit = qubit.name
         else:
             pulse = platform.create_qubit_flux_pulse(
                 qubit, start=0, duration=sequence.duration
             )
         qf_pulses[qubit] = pulse
-        sequence.add(pulse)
+        if crosstalk:
+            sequences[i].add(pulse)
+        else:
+            sequence.add(pulse)
 
-    # FIXME: This is a patch to fix couplers/qubits
-    if isinstance(qubits[0], Coupler):
-        pulses = [qf_pulses[qubit.name] for qubit in qubits]
+    if crosstalk:
+        sweepers = [
+            Sweeper(
+                Parameter.amplitude,
+                delta_bias_flux_range,
+                pulses=[qf_pulses[qubit]],
+                type=SweeperType.ABSOLUTE,
+            )
+            for qubit in qubits
+        ]
+        return delta_bias_flux_range, sweepers, sequences
     else:
-        pulses = [qf_pulses[qubit] for qubit in qubits]
-
-    sweepers = [
-        Sweeper(
-            Parameter.amplitude,
-            delta_bias_flux_range,
-            pulses=pulses,
-            type=SweeperType.ABSOLUTE,
-        )
-    ]
-    return delta_bias_flux_range, sweepers
+        sweepers = [
+            Sweeper(
+                Parameter.amplitude,
+                delta_bias_flux_range,
+                pulses=[qf_pulses[qubit] for qubit in qubits],
+                type=SweeperType.ABSOLUTE,
+            )
+        ]
+        return delta_bias_flux_range, sweepers, [sequence]
 
 
 def _acquisition(
@@ -228,12 +239,13 @@ def _acquisition(
     ro_pulses = {}
     qubit_frequency = {}
     bare_resonator_frequency = {}
+    offset = {}
     for qubit in targets:
         qubit_frequency[qubit] = platform.qubits[qubit].drive_frequency
         bare_resonator_frequency[qubit] = platform.qubits[
             qubit
         ].bare_resonator_frequency
-
+        offset[qubit] = platform.qubits[qubit].sweetspot
         ro_pulses[qubit] = platform.create_qubit_readout_pulse(qubit, start=0)
         sequence.add(ro_pulses[qubit])
 
@@ -248,9 +260,10 @@ def _acquisition(
         type=SweeperType.OFFSET,
     )
     if params.flux_pulses:
-        delta_bias_flux_range, sweepers = create_flux_pulse_sweepers(
+        delta_bias_flux_range, sweepers, sequences = create_flux_pulse_sweepers(
             params, platform, targets, sequence
         )
+        sequence = sequences[0]
     else:
         delta_bias_flux_range = np.arange(
             -params.bias_width / 2, params.bias_width / 2, params.bias_step
@@ -268,6 +281,7 @@ def _acquisition(
         resonator_type=platform.resonator_type,
         flux_pulses=params.flux_pulses,
         qubit_frequency=qubit_frequency,
+        offset=offset,
         bare_resonator_frequency=bare_resonator_frequency,
     )
 
@@ -321,7 +335,7 @@ def _fit(data: ResonatorFluxData) -> ResonatorFluxResults:
 
         try:
             popt = curve_fit(
-                utils.transmon_readout_frequency,
+                utils.transmon_readout_frequency_diagonal,
                 biases,
                 frequencies * HZ_TO_GHZ,
                 bounds=utils.resonator_flux_dependence_fit_bounds(
@@ -336,7 +350,7 @@ def _fit(data: ResonatorFluxData) -> ResonatorFluxResults:
             # frequency corresponds to transmon readout frequency
             # at the sweetspot popt[3]
             frequency[qubit] = (
-                utils.transmon_readout_frequency(popt[3], *popt) * GHZ_TO_HZ
+                utils.transmon_readout_frequency_diagonal(popt[3], *popt) * GHZ_TO_HZ
             )
             sweetspot[qubit] = popt[3]
             asymmetry[qubit] = popt[1]
@@ -367,7 +381,7 @@ def _fit(data: ResonatorFluxData) -> ResonatorFluxResults:
 def _plot(data: ResonatorFluxData, fit: ResonatorFluxResults, target: QubitId):
     """Plotting function for ResonatorFlux Experiment."""
     figures = utils.flux_dependence_plot(
-        data, fit, target, utils.transmon_readout_frequency
+        data, fit, target, utils.transmon_readout_frequency_diagonal
     )
     if data.flux_pulses:
         bias_flux_unit = "a.u."
@@ -384,7 +398,7 @@ def _plot(data: ResonatorFluxData, fit: ResonatorFluxResults, target: QubitId):
                     "Qubit Frequency at Sweetspot [Hz]",
                     "Asymmetry d",
                     "Coupling g",
-                    "V_ii [V]",
+                    "Flux dependence",
                 ],
                 [
                     np.round(fit.sweetspot[target], 4),
diff --git a/src/qibocal/protocols/characterization/flux_dependence/utils.py b/src/qibocal/protocols/characterization/flux_dependence/utils.py
index 3127eaee9..3b0e32fb7 100644
--- a/src/qibocal/protocols/characterization/flux_dependence/utils.py
+++ b/src/qibocal/protocols/characterization/flux_dependence/utils.py
@@ -48,7 +48,11 @@ def flux_dependence_plot(data, fit, qubit, fit_function=None):
     fig.add_trace(
         go.Heatmap(
             x=qubit_data.freq * HZ_TO_GHZ,
-            y=qubit_data.bias,
+            y=(
+                qubit_data.bias - data.offset[qubit]
+                if data.flux_pulses
+                else qubit_data.bias
+            ),
             z=qubit_data.signal,
             colorbar_x=0.46,
         ),
@@ -67,10 +71,10 @@ def flux_dependence_plot(data, fit, qubit, fit_function=None):
         fig.add_trace(
             go.Scatter(
                 x=fit_function(bias, *params),
-                y=bias,
+                y=bias - data.offset[qubit] if data.flux_pulses else bias,
                 showlegend=True,
                 name="Fit",
-                marker=dict(color="black"),
+                marker=dict(color="green"),
             ),
             row=1,
             col=1,
@@ -83,7 +87,7 @@ def flux_dependence_plot(data, fit, qubit, fit_function=None):
     )
 
     if data.flux_pulses:
-        fig.update_yaxes(title_text="Flux [a.u.]", row=1, col=1)
+        fig.update_yaxes(title_text="Flux Amplitude [a.u.]", row=1, col=1)
     else:
         fig.update_yaxes(title_text="Bias [V]", row=1, col=1)
 
@@ -115,16 +119,14 @@ def flux_dependence_plot(data, fit, qubit, fit_function=None):
     return figures
 
 
-def flux_crosstalk_plot(data, qubit):
+def flux_crosstalk_plot(data, qubit, fit, fit_function):
     figures = []
     fitting_report = ""
-
     all_qubit_data = {
         index: data_qubit
         for index, data_qubit in data.data.items()
         if index[0] == qubit
     }
-
     fig = make_subplots(
         rows=1,
         cols=len(all_qubit_data),
@@ -134,19 +136,67 @@ def flux_crosstalk_plot(data, qubit):
     )
     for col, (flux_qubit, qubit_data) in enumerate(all_qubit_data.items()):
         frequencies = qubit_data.freq * HZ_TO_GHZ
-        signal = qubit_data.signal
-        if data.resonator_type == "2D":
-            signal = -signal
-
         fig.add_trace(
             go.Heatmap(
                 x=frequencies,
-                y=qubit_data.bias,
+                y=(
+                    qubit_data.bias - data.voltage[qubit]
+                    if data.flux_pulses
+                    else qubit_data.bias
+                ),
                 z=qubit_data.signal,
+                showscale=False,
             ),
             row=1,
             col=col + 1,
         )
+        if fit is not None:
+
+            if flux_qubit[1] != qubit:
+                fig.add_trace(
+                    go.Scatter(
+                        x=fit_function(
+                            xj=qubit_data.bias, **fit.fitted_parameters[flux_qubit]
+                        )
+                        * HZ_TO_GHZ,
+                        y=(
+                            qubit_data.bias - data.voltage[qubit]
+                            if data.flux_pulses
+                            else qubit_data.bias
+                        ),
+                        showlegend=not any(
+                            isinstance(trace, go.Scatter) for trace in fig.data
+                        ),
+                        legendgroup="Fit",
+                        name="Fit",
+                        marker=dict(color="green"),
+                    ),
+                    row=1,
+                    col=col + 1,
+                )
+            else:
+                diagonal_params = fit.fitted_parameters[qubit, qubit]
+                fig.add_trace(
+                    go.Scatter(
+                        x=globals().get(fit_function.__name__ + "_diagonal")(
+                            qubit_data.bias,
+                            *diagonal_params,
+                        ),
+                        y=(
+                            qubit_data.bias - data.voltage[flux_qubit[1]]
+                            if data.flux_pulses
+                            else qubit_data.bias
+                        ),
+                        showlegend=not any(
+                            isinstance(trace, go.Scatter) for trace in fig.data
+                        ),
+                        legendgroup="Fit",
+                        name="Fit",
+                        marker=dict(color="green"),
+                    ),
+                    row=1,
+                    col=col + 1,
+                )
 
         fig.update_xaxes(
             title_text="Frequency [GHz]",
@@ -156,7 +206,9 @@ def flux_crosstalk_plot(data, qubit):
 
         if data.flux_pulses:
             fig.update_yaxes(
-                title_text=f"Qubit {flux_qubit[1]}: Flux [a.u.]", row=1, col=col + 1
+                title_text=f"Qubit {flux_qubit[1]}: Flux Amplitude [a.u.]",
+                row=1,
+                col=col + 1,
             )
         else:
             fig.update_yaxes(
@@ -164,35 +216,46 @@ def flux_crosstalk_plot(data, qubit):
             )
 
     fig.update_layout(xaxis1=dict(range=[np.min(frequencies), np.max(frequencies)]))
-    fig.update_traces(showscale=False)  # disable colorbar
+    fig.update_layout(xaxis2=dict(range=[np.min(frequencies), np.max(frequencies)]))
+    fig.update_layout(xaxis3=dict(range=[np.min(frequencies), np.max(frequencies)]))
     fig.update_layout(
-        showlegend=False,
+        showlegend=True,
     )
-
     figures.append(fig)
 
     return figures, fitting_report
 
 
-def G_f_d(x, offset, d, element):
+def G_f_d(xi, xj, sweetspot, d, matrix_element, crosstalk_element):
     """Auxiliary function to calculate qubit frequency as a function of bias.
 
-    It also determines the flux dependence of :math:`E_J`,:math:`E_J(\\phi)=E_J(0)G_f_d^2`.
+    It also determines the flux dependence of :math:`E_J`,:math:`E_J(\\phi)=E_J(0)G_f_d`.
     For more details see: https://arxiv.org/pdf/cond-mat/0703002.pdf
 
     Args:
-        offset (float): bias offset.
-        matrix_element(float): constant to convert flux (:math:`\\phi_0`) to bias (:math:`v_0`). Typically denoted as :math:`\\Xi`. :math:`v_0 = \\Xi \\phi_0`.
+        xi (float): bias of target qubit
+        xj (float): bias of neighbor qubit
+        sweetspot (float): sweetspot [V].
+        matrix_element(float): diagonal crosstalk matrix element
+        crosstalk_element(float): off-diagonal crosstalk matrix element
         d (float): asymmetry between the two junctions of the transmon.
                    Typically denoted as :math:`d`. :math:`d = (E_J^1 - E_J^2) / (E_J^1 + E_J^2)`.
 
     Returns:
         (float)
     """
-    return (d**2 + (1 - d**2) * np.cos(np.pi * (x - offset) * element) ** 2) ** 0.25
+    return (
+        d**2
+        + (1 - d**2)
+        * np.cos(
+            np.pi
+            * ((xi - sweetspot) * matrix_element + (xj - sweetspot) * crosstalk_element)
+        )
+        ** 2
+    ) ** 0.25
 
 
-def transmon_frequency(x, w_max, d, element, offset):
+def transmon_frequency(xi, xj, w_max, d, matrix_element, sweetspot, crosstalk_element):
     r"""Approximation to transmon frequency.
 
     The formula holds in the transmon regime Ej / Ec >> 1.
@@ -200,39 +263,133 @@ def transmon_frequency(x, w_max, d, element, offset):
     See  https://arxiv.org/pdf/cond-mat/0703002.pdf for the complete formula.
 
     Args:
-         x (float): bias value
-         w_max (float): maximum frequency  :math:`w_{max} = \sqrt{8 E_j E_c}
-         d (float):  d (float): asymmetry between the two junctions of the transmon.
-         element (float): matrix element
-         offset (float): bias corresponding to zero flux (sweetspot).
+        xi (float): bias of target qubit
+        xj (float): bias of neighbor qubit
+        w_max (float): maximum frequency  :math:`w_{max} = \sqrt{8 E_j E_c}
+        sweetspot (float): sweetspot [V].
+        matrix_element(float): diagonal crosstalk matrix element
+        crosstalk_element(float): off-diagonal crosstalk matrix element
+        d (float): asymmetry between the two junctions of the transmon.
+                   Typically denoted as :math:`d`. :math:`d = (E_J^1 - E_J^2) / (E_J^1 + E_J^2)`.
 
      Returns:
          (float): qubit frequency as a function of bias.
     """
-    return w_max * G_f_d(x, offset=offset, d=d, element=element)
+    return w_max * G_f_d(
+        xi,
+        xj,
+        sweetspot=sweetspot,
+        d=d,
+        matrix_element=matrix_element,
+        crosstalk_element=crosstalk_element,
+    )
+
+
+def transmon_frequency_diagonal(x, w_max, d, matrix_element, sweetspot):
+    r"""Approximation to transmon frequency without crosstalk.
+
+    The formula holds in the transmon regime Ej / Ec >> 1.
 
+    See  https://arxiv.org/pdf/cond-mat/0703002.pdf for the complete formula.
+
+    Args:
+        x (float): bias of target qubit
+        w_max (float): maximum frequency  :math:`w_{max} = \sqrt{8 E_j E_c}
+        sweetspot (float): sweetspot [V].
+        matrix_element(float): diagonal crosstalk matrix element
+        d (float): asymmetry between the two junctions of the transmon.
+                   Typically denoted as :math:`d`. :math:`d = (E_J^1 - E_J^2) / (E_J^1 + E_J^2)`.
 
-def transmon_readout_frequency(x, w_max, d, element, offset, resonator_freq, g):
+     Returns:
+         (float): qubit frequency as a function of bias.
+    """
+    return transmon_frequency(
+        xi=x,
+        xj=0,
+        w_max=w_max,
+        d=d,
+        matrix_element=matrix_element,
+        sweetspot=sweetspot,
+        crosstalk_element=0,
+    )
+
+
+def transmon_readout_frequency(
+    xi, xj, w_max, d, matrix_element, crosstalk_element, sweetspot, resonator_freq, g
+):
     r"""Approximation to flux dependent resonator frequency.
 
     The formula holds in the transmon regime Ej / Ec >> 1.
 
     See  https://arxiv.org/pdf/cond-mat/0703002.pdf for the complete formula.
 
+    Args:
+         xi (float): bias of target qubit
+         xj (float): bias of neighbor qubit
+         w_max (float): maximum frequency  :math:`w_{max} = \sqrt{8 E_j E_c}
+         sweetspot (float): sweetspot [V].
+         matrix_element(float): diagonal crosstalk matrix element
+         crosstalk_element(float): off-diagonal crosstalk matrix element
+         d (float): asymmetry between the two junctions of the transmon.
+                    Typically denoted as :math:`d`. :math:`d = (E_J^1 - E_J^2) / (E_J^1 + E_J^2)`.
+         resonator_freq (float): bare resonator frequency [GHz]
+         g (float): readout coupling.
+
+     Returns:
+         (float): resonator frequency as a function of bias.
+    """
+    return resonator_freq + g**2 * G_f_d(
+        xi=xi,
+        xj=xj,
+        sweetspot=sweetspot,
+        d=d,
+        matrix_element=matrix_element,
+        crosstalk_element=crosstalk_element,
+    ) / (
+        resonator_freq
+        - transmon_frequency(
+            xi=xi,
+            xj=xj,
+            w_max=w_max,
+            d=d,
+            matrix_element=matrix_element,
+            sweetspot=sweetspot,
+            crosstalk_element=crosstalk_element,
+        )
+    )
+
+
+def transmon_readout_frequency_diagonal(
+    x, w_max, d, matrix_element, sweetspot, resonator_freq, g
+):
+    r"""Approximation to flux dependent resonator frequency without crosstalk.
+
+    The formula holds in the transmon regime Ej / Ec >> 1.
+
+    See  https://arxiv.org/pdf/cond-mat/0703002.pdf for the complete formula.
+
     Args:
          x (float): bias value
          w_max (float): maximum frequency  :math:`w_{max} = \sqrt{8 E_j E_c}
          d (float):  d (float): asymmetry between the two junctions of the transmon.
-         element (float): matrix element
-         offset (float): bias corresponding to zero flux (sweetspot).
+         matrix element (float): diagonal crosstalk matrix element
+         phase (float): bias corresponding to zero flux (sweetspot).
          resonator_freq (float): bare resonator frequency [GHz]
          g (float): readout coupling.
 
      Returns:
          (float): resonator frequency as a function of bias.
     """
-    return resonator_freq + g**2 * G_f_d(x, offset, d, element) / (
-        resonator_freq - transmon_frequency(x, w_max, d, element, offset)
+    return transmon_readout_frequency(
+        xi=x,
+        xj=0,
+        w_max=w_max,
+        sweetspot=sweetspot,
+        d=d,
+        matrix_element=matrix_element,
+        crosstalk_element=0,
+        g=g,
+        resonator_freq=resonator_freq,
     )
 
 
diff --git a/src/qibocal/update.py b/src/qibocal/update.py
index 584bf2c76..7d197d46e 100644
--- a/src/qibocal/update.py
+++ b/src/qibocal/update.py
@@ -65,6 +65,13 @@ def drive_duration(duration: Union[int, tuple], platform: Platform, qubit: Qubit
     platform.qubits[qubit].native_gates.RX.duration = int(duration)
 
 
+def crosstalk_matrix(
+    matrix_element: float, platform: Platform, qubit: QubitId, flux_qubit: QubitId
+):
+    """Update crosstalk_matrix element."""
+    platform.qubits[qubit].crosstalk_matrix[flux_qubit] = float(matrix_element)
+
+
 def iq_angle(angle: float, platform: Platform, qubit: QubitId):
     """Update iq angle value in platform for specific qubit."""
     platform.qubits[qubit].iq_angle = float(angle)