Evaluate the sweetspot closest to the center of bias interval

doc/source/protocols/flux/single.rst

@@ -16,9 +16,8 @@ of GHz, which leads to several applications including quantum logical gates.
 
 The transmon frequency as a function of the external flux can be expressed as :cite:p:`Barrett_2023`
 
-.. math::
-
-    f_q(\Phi) = \Bigg( f_q^{\text{max}} + \frac{E_C}{h} \Bigg)  \sqrt[4]{d^2 + (1-d^2)\cos^2\Big( \pi \frac{\Phi}{\Phi_0}\Big)} - \frac{E_C}{h} \,
+.. math:: f_q(\Phi) = \Bigg( f_q^{\text{max}} + \frac{E_C}{h} \Bigg)  \sqrt[4]{d^2 + (1-d^2)\cos^2\Big( \pi \frac{\Phi}{\Phi_0}\Big)} - \frac{E_C}{h} \,
+    :label: transmon
 
 where :math:`f_{\text{max}} = ( \sqrt{8 E_C E_J} - E_C) / h` is the maximum qubit frequency,
 :math:`d` is the junctions asymmetry, :math:`E_C` is the charging energy,
@@ -76,6 +75,14 @@ The expected output is the following:
 From the acquired data this protocol estimates the flux insensitive point "sweetspot",
 which corresponds to the flux value where the frequency is maximed, as well as the drive frequency
 and the diagonal crosstalk coefficient :math:`V_{ii}`.
+
+.. note::
+
+    From the cosinusoidal term in the transmon equation :math:numref:`transmon`, it is clear that the
+    sweetspot is not unique.
+    In this protocol, Qibocal returns the sweetspot that is closest to the bias
+    that is in the middle of the swept interval.
+
 Both the sweetspot and the :math:`C_{ii}` can be understood by writing the full expression for
 the flux felt by qubit :math:`i` :cite:p:`Barrett_2023`:
 

src/qibocal/protocols/flux_dependence/qubit_flux_dependence.py

@@ -165,9 +165,14 @@ def _acquisition(
 
 def _fit(data: QubitFluxData) -> QubitFluxResults:
     """
-    Post-processing for QubitFlux Experiment. See arxiv:0703002
-    Fit frequency as a function of current for the flux qubit spectroscopy
-    data (QubitFluxData): data object with information on the feature response at each current point.
+    Post-processing for QubitFlux Experiment. See `arXiv:0703002 <https://arxiv.org/abs/cond-mat/0703002>`_.
+    Fit frequency as a function of current for the flux qubit spectroscopy data.
+    All possible sweetspots :math:`x` are evaluated by the function
+    :math:`x p_1 + p_2 = k`, for integers :math:`k`, where :math:`p_1` and :math:`p_2`
+    are respectively the normalization and the offset, as defined in
+    :mod:`qibocal.protocols.flux_dependence.utils.transmon_frequency`.
+    The code returns the sweetspot that is closest to the bias
+    in the middle of the swept interval.
     """
 
     qubits = data.qubits
@@ -178,12 +183,15 @@ def _fit(data: QubitFluxData) -> QubitFluxResults:
 
     for qubit in qubits:
         qubit_data = data[qubit]
-        biases = qubit_data.bias
-        frequencies = qubit_data.freq
+        interval_biases = qubit_data.bias
+        interval_frequencies = qubit_data.freq
         signal = qubit_data.signal
 
         frequencies, biases = extract_feature(
-            frequencies, biases, signal, "max" if data.resonator_type == "2D" else "min"
+            interval_frequencies,
+            interval_biases,
+            signal,
+            "max" if data.resonator_type == "2D" else "min",
         )
 
         def fit_function(x, w_max, normalization, offset):
@@ -218,10 +226,10 @@ def fit_function(x, w_max, normalization, offset):
                 "charging_energy": data.charging_energy[qubit] * HZ_TO_GHZ,
             }
             frequency[qubit] = popt[0] * GHZ_TO_HZ
-            # solution to x*popt[1] + popt[2] = k
-            # such that x is close to 0
-            # to avoid errors due to periodicity
-            sweetspot[qubit] = (np.round(popt[2]) - popt[2]) / popt[1]
+            middle_bias = (np.max(interval_biases) + np.min(interval_biases)) / 2
+            sweetspot[qubit] = (
+                np.round(popt[1] * middle_bias + popt[2]) - popt[2]
+            ) / popt[1]
             matrix_element[qubit] = popt[1]
         except ValueError as e:
             log.error(

src/qibocal/protocols/flux_dependence/utils.py

@@ -205,11 +205,10 @@ def G_f_d(xi, xj, offset, d, crosstalk_element, normalization):
         xi (float): bias of target qubit
         xj (float): bias of neighbor qubit
         offset (float): phase_offset [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)`.
-        normalization (float): Normalize diagonal element to 1
+        crosstalk_element(float): off-diagonal crosstalk matrix element
+        normalization(float): diagonal crosstalk matrix element
     Returns:
         (float)
     """
@@ -237,11 +236,11 @@ def transmon_frequency(
         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)`.
+        normalization(float): diagonal crosstalk matrix element
+        offset (float): phase_offset [V].
+        crosstalk_element(float): off-diagonal crosstalk matrix element
         charging_energy (float): Ec / h (GHz)
 
      Returns:
@@ -279,13 +278,14 @@ def transmon_readout_frequency(
          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)`.
+         normalization(float): diagonal crosstalk matrix element
+         offset (float): phase_offset [V].
+         crosstalk_element(float): off-diagonal crosstalk matrix element
          resonator_freq (float): bare resonator frequency [GHz]
          g (float): readout coupling.
+         charging_energy (float): Ec / h (GHz)
 
      Returns:
          (float): resonator frequency as a function of bias.