diff --git a/docs/index.rst b/docs/index.rst
index f49ec2c..ae01ae7 100755
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -19,15 +19,19 @@ functionals or data formats.
 Contributors
 ============
 
-- **Monika Marek** (2023)
+- **Monika Marek** (2023) - `mm1707 <https://github.com/mm1707>`_
 
   Developed module :mod:`.tools`.
 
-- **Aleksandra Bochenek** (2022)
+- **Aleksandra Bochenek** (2022) - `pirrx <https://github.com/pirrx>`_
 
   Developed :mod:`.bsk` by implementing and testing equations for the Brussels-Montreal family
   of density functionals.
 
+- **Adarsh Karekkat** (2025) - `O-slyfox-O <https://github.com/O-slyfox-O>`_
+
+  Added pairing interpolation schemes to :mod:`.bsk` for newer Bsk models including BskG4.
+
 Acknowledgments
 ===============
 Funding
diff --git a/libnest/bsk.py b/libnest/bsk.py
index 542e93c..219669a 100644
--- a/libnest/bsk.py
+++ b/libnest/bsk.py
@@ -339,7 +339,7 @@ def proton_ref_pairing_field(rho_n, rho_p):
     rho, eta = rhoEta(rho_n, rho_p) #eta = rho_n - rho_p
     rho = np.asarray(rho + DENSEPSILON, dtype=float)
     return (symmetric_pairing_field(rho_n, rho_p)*(1-np.abs(eta/rho))
-        -neutron_pairing_field(rho_n)*rho_p/rho*eta/rho)
+        -neutron_pairing_field(rho_n)*rho_p/rho*eta/rho)                #Adarsh: should it be rho_p here in neutron_pairing_field()?
 
 
 # ================================
@@ -351,23 +351,261 @@ def proton_ref_pairing_field(rho_n, rho_p):
 # (neutron_pairing_field). Add two more schemes here, following:
 # https://link.springer.com/article/10.1140/epja/s10050-025-01503-x#citeas
 
-def ref_pairing_field_eq2(rho_n, rho_p):
-    """TODO(Adarsh): implement the interpolation scheme from Eq. (2) of
+def neutron_ref_pairing_field_eq2(rho_n, rho_p):
+    r"""
+    Returns the reference pairing field for neutrons in uniform matter
+    using the interpolation scheme from Eq. (2) of
     https://link.springer.com/article/10.1140/epja/s10050-025-01503-x#citeas
+
+    This interpolates between symmetric matter and neutron matter,
+    with asymmetry dependence governed by :math:`\delta`.
+    
+    Used in all Bsk parameterizations since Bsk17. Currently in use in WBSK.
+
+    .. math::
+
+        \Delta_n(\rho_n,\rho_p) =
+        \Delta_{\mathrm{SM}}(\rho)
+        \left(1 - |\delta|\right)
+        + \delta \frac{\rho_p}{\rho}
+        \Delta_{\mathrm{NeuM}}(\rho_n)
+
+    where :math:`\rho = \rho_n + \rho_p` and
+    :math:`\delta = (\rho_n - \rho_p)/\rho`.
+
+    Args:
+        rho_n (float): neutron density :math:`\rho_n` [fm :sup:`-3`]; sum of both spin components
+        rho_p (float): proton density :math:`\rho_p` [fm :sup:`-3`]; sum of both spin components
+
+    Returns:
+        float: pairing field for neutrons :math:`\Delta_n` [MeV]
+
+    See also:
+        :func:`.proton_ref_pairing_field_eq2`
+        :func:`.symmetric_pairing_field`
+        :func:`.neutron_pairing_field`
     """
-    raise NotImplementedError
+    
+    rho_n = np.asarray(rho_n, dtype=float)
+    rho_p = np.asarray(rho_p, dtype=float)
+    rho, eta = rhoEta(rho_n, rho_p) #eta = rho_n - rho_p
+    rho = np.asarray(rho + DENSEPSILON, dtype=float)
+    delta = eta/rho
+    pairing = ((1-np.abs(delta))*symmetric_pairing_field(rho_n,rho_p)) + (delta* (rho_n/rho) * neutron_pairing_field(rho_n))
+    return(pairing)
 
-def ref_pairing_field_eq3(rho_n, rho_p):
-    """TODO(Adarsh): implement the interpolation scheme from Eq. (3) of
+def proton_ref_pairing_field_eq2(rho_n, rho_p):
+    
+    r"""
+    Returns the reference pairing field for protons in uniform matter
+    using the interpolation scheme from Eq. (2) of
     https://link.springer.com/article/10.1140/epja/s10050-025-01503-x#citeas
+
+    This interpolates between symmetric matter and neutron matter,
+    with asymmetry dependence governed by :math:`\delta`.
+    
+    Used in all Bsk parameterizations since Bsk17. Currently in use in WBSK.
+
+    .. math::
+
+        \Delta_p(\rho_n,\rho_p) =
+        \Delta_{\mathrm{SM}}(\rho)
+        \left(1 - |\delta|\right)
+        - \delta \frac{\rho_p}{\rho}
+        \Delta_{\mathrm{NeuM}}(\rho_p)
+
+    where :math:`\rho = \rho_n + \rho_p` and
+    :math:`\delta = (\rho_n - \rho_p)/\rho`.
+
+    Args:
+        rho_n (float): neutron density :math:`\rho_n` [fm :sup:`-3`]
+        rho_p (float): proton density :math:`\rho_p` [fm :sup:`-3`]
+
+    Returns:
+        float: pairing field for protons :math:`\Delta_p` [MeV]
+
+    See also:
+        :func:`.neutron_ref_pairing_field_eq2`
+        :func:`.symmetric_pairing_field`
+        :func:`.neutron_pairing_field`
     """
-    raise NotImplementedError
+    rho_n = np.asarray(rho_n, dtype=float)
+    rho_p = np.asarray(rho_p, dtype=float)
+    rho, eta = rhoEta(rho_n, rho_p) #eta = rho_n - rho_p
+    rho = np.asarray(rho + DENSEPSILON, dtype=float)
+    delta = eta/rho
+    pairing = ((1-np.abs(delta))*symmetric_pairing_field(rho_n,rho_p)) - (delta* (rho_p/rho) * neutron_pairing_field(rho_p))
+    return(pairing)
 
-def ref_pairing_field_eq6(rho_n, rho_p):
-    """TODO(Adarsh): implement the interpolation scheme from Eq. (6) of
+def neutron_ref_pairing_field_eq3(rho_n, rho_p):
+    
+    r"""
+    Returns the reference pairing field for neutrons in uniform matter
+    using the interpolation scheme from Eq. (3) of
     https://link.springer.com/article/10.1140/epja/s10050-025-01503-x#citeas
+
+    This interpolates between symmetric matter and neutron matter,
+    with asymmetry dependence governed by :math:`\delta`.
+    
+    Used in BskG3. 
+
+    .. math::
+
+        \Delta_n(\rho_n,\rho_p) =
+        \Delta_{\mathrm{SM}}(\rho)
+        \left(1 - |\delta|\right)
+        + |\delta| \,
+        \Delta_{\mathrm{NeuM}}(\rho_n)
+
+    where :math:`\delta = (\rho_n - \rho_p)/\rho`.
+
+    Args:
+        rho_n (float): neutron density :math:`\rho_n` [fm :sup:`-3`]
+        rho_p (float): proton density :math:`\rho_p` [fm :sup:`-3`]
+
+    Returns:
+        float: pairing field for neutrons :math:`\Delta_n` [MeV]
+
+    See also:
+        :func:`.proton_ref_pairing_field_eq3`
+        :func:`.symmetric_pairing_field`
+        :func:`.neutron_pairing_field`
     """
-    raise NotImplementedError
+    rho_n = np.asarray(rho_n, dtype=float)
+    rho_p = np.asarray(rho_p, dtype=float)
+    rho, eta = rhoEta(rho_n, rho_p) #eta = rho_n - rho_p
+    rho = np.asarray(rho + DENSEPSILON, dtype=float)
+    delta = eta/rho
+    pairing = ((1-np.abs(delta))*symmetric_pairing_field(rho_n,rho_p)) + np.abs(delta) * neutron_pairing_field(rho_n)
+    return(pairing)
+
+def proton_ref_pairing_field_eq3(rho_n, rho_p):
+
+    r"""
+    Returns the reference pairing field for protons in uniform matter
+    using the interpolation scheme from Eq. (3) of
+    https://link.springer.com/article/10.1140/epja/s10050-025-01503-x#citeas
+
+    This interpolates between symmetric matter and neutron matter,
+    with asymmetry dependence governed by :math:`\delta`.
+    
+    Used in BskG3.
+
+    .. math::
+
+        \Delta_p(\rho_n,\rho_p) =
+        \Delta_{\mathrm{SM}}(\rho)
+        \left(1 - |\delta|\right)
+        + |\delta| \,
+        \Delta_{\mathrm{NeuM}}(\rho_p)
+
+    where :math:`\delta = (\rho_n - \rho_p)/\rho`.
+
+    Args:
+        rho_n (float): neutron density :math:`\rho_n` [fm :sup:`-3`]
+        rho_p (float): proton density :math:`\rho_p` [fm :sup:`-3`]
+
+    Returns:
+        float: pairing field for protons :math:`\Delta_p` [MeV]
+
+    See also:
+        :func:`.neutron_ref_pairing_field_eq3`
+        :func:`.symmetric_pairing_field`
+        :func:`.neutron_pairing_field`
+    """
+    rho_n = np.asarray(rho_n, dtype=float)
+    rho_p = np.asarray(rho_p, dtype=float)
+    rho, eta = rhoEta(rho_n, rho_p) #eta = rho_n - rho_p
+    rho = np.asarray(rho + DENSEPSILON, dtype=float)
+    delta = eta/rho
+    pairing = ((1-np.abs(delta))*symmetric_pairing_field(rho_n,rho_p)) + np.abs(delta) * neutron_pairing_field(rho_p)
+    return(pairing)
+
+def neutron_ref_pairing_field_eq6(rho_n, rho_p):
+
+    r"""
+    Returns the reference pairing field for neutrons in uniform matter
+    using the interpolation scheme from Eq. (6) of
+    https://link.springer.com/article/10.1140/epja/s10050-025-01503-x#citeas
+
+    This interpolates between symmetric matter and neutron matter,
+    with asymmetry dependence governed by :math:`\delta`.
+    
+    Used in BskG4.
+
+    .. math::
+
+        \Delta_n(\rho_n,\rho_p) =
+        \Delta_{\mathrm{NeuM}}(\rho_n)
+        \left[
+        \frac{\Delta_{\mathrm{SM}}(\rho)}
+             {\Delta_{\mathrm{NeuM}}(\rho/2)}
+        \right]^{(1 - \delta)}
+
+    where :math:`\delta = (\rho_n - \rho_p)/\rho`.
+
+    Args:
+        rho_n (float): neutron density :math:`\rho_n` [fm :sup:`-3`]
+        rho_p (float): proton density :math:`\rho_p` [fm :sup:`-3`]
+
+    Returns:
+        float: pairing field for neutrons :math:`\Delta_n` [MeV]
+
+    See also:
+        :func:`.proton_ref_pairing_field_eq6`
+        :func:`.symmetric_pairing_field`
+        :func:`.neutron_pairing_field`
+    """
+    rho_n = np.asarray(rho_n, dtype=float)
+    rho_p = np.asarray(rho_p, dtype=float)
+    rho, eta = rhoEta(rho_n, rho_p) #eta = rho_n - rho_p
+    rho = np.asarray(rho + DENSEPSILON, dtype=float)
+    delta = eta/rho
+    pairing = neutron_pairing_field(rho_n)* ((symmetric_pairing_field(rho_n,rho_p)/neutron_pairing_field(rho/2))**(1-delta))
+    return(pairing)
+
+def proton_ref_pairing_field_eq6(rho_n, rho_p):
+
+    r"""
+    Returns the reference pairing field for protons in uniform matter
+    using the interpolation scheme from Eq. (6) of
+    https://link.springer.com/article/10.1140/epja/s10050-025-01503-x#citeas
+
+    This interpolates between symmetric matter and neutron matter,
+    with asymmetry dependence governed by :math:`\delta`.
+    
+    Used in BskG4.
+
+    .. math::
+
+        \Delta_p(\rho_n,\rho_p) =
+        \Delta_{\mathrm{NeuM}}(\rho_p)
+        \left[
+        \frac{\Delta_{\mathrm{SM}}(\rho)}
+             {\Delta_{\mathrm{NeuM}}(\rho/2)}
+        \right]^{(1 + \delta)}
+
+    where :math:`\delta = (\rho_n - \rho_p)/\rho`.
+
+    Args:
+        rho_n (float): neutron density :math:`\rho_n` [fm :sup:`-3`]
+        rho_p (float): proton density :math:`\rho_p` [fm :sup:`-3`]
+
+    Returns:
+        float: pairing field for protons :math:`\Delta_p` [MeV]
+
+    See also:
+        :func:`.neutron_ref_pairing_field_eq6`
+        :func:`.symmetric_pairing_field`
+        :func:`.neutron_pairing_field`
+    """
+    rho_n = np.asarray(rho_n, dtype=float)
+    rho_p = np.asarray(rho_p, dtype=float)
+    rho, eta = rhoEta(rho_n, rho_p) #eta = rho_n - rho_p
+    rho = np.asarray(rho + DENSEPSILON, dtype=float)
+    delta = eta/rho
+    pairing = neutron_pairing_field(rho_p)* ((symmetric_pairing_field(rho_n,rho_p)/neutron_pairing_field(rho/2))**(1+delta))
+    return(pairing)
 
 
 # ================================
diff --git a/tests/test_bsk.py b/tests/test_bsk.py
index 2f3879f..bf2d4ec 100644
--- a/tests/test_bsk.py
+++ b/tests/test_bsk.py
@@ -212,36 +212,122 @@ class TestAlternativePairingSchemes(unittest.TestCase):
     """Tests for the alternative delta_n/delta_p interpolation schemes
     (bsk.ref_pairing_field_eq2/eq3/eq6), based on Eq. (2), (3) and (6) of
     https://link.springer.com/article/10.1140/epja/s10050-025-01503-x#citeas
-
-    TODO(Adarsh): each scheme is currently a stub that raises
-    NotImplementedError (see bsk.py). As you implement a scheme, remove its
-    @unittest.skip decorator and fill in real assertions - e.g. following the
-    pattern in TestReferencePairingFields above:
-      - eta = 0 (rho_n == rho_p) should reduce to symmetric_pairing_field
-      - rho_p = 0 should reduce to neutron_pairing_field
-      - finite output for both scalar and array input
     """
 
-    @unittest.skip("TODO(Adarsh): implement ref_pairing_field_eq2 (Eq. 2)")
-    def test_ref_pairing_field_eq2(self):
+    def test_neutron_ref_pairing_field_eq2(self):
+        # Symmetric limit
         rho_n, rho_p = 0.05, 0.05
-        ref = bsk.ref_pairing_field_eq2(rho_n, rho_p)
+        ref = bsk.neutron_ref_pairing_field_eq2(rho_n, rho_p)
         sym = bsk.symmetric_pairing_field(rho_n, rho_p)
         self.assertAlmostEqual(float(ref), float(sym), places=6)
-
-    @unittest.skip("TODO(Adarsh): implement ref_pairing_field_eq3 (Eq. 3)")
-    def test_ref_pairing_field_eq3(self):
+        # Pure neutron matter limit
+        rho_n, rho_p = 0.05, 0.00
+        ref  = bsk.neutron_ref_pairing_field_eq2(rho_n, rho_p)
+        neum = bsk.neutron_pairing_field(rho_n)
+        self.assertAlmostEqual(float(ref), float(neum), places=6)
+        # Pure proton matter limit
+        rho_n, rho_p = 0.00, 0.05
+        ref   = bsk.neutron_ref_pairing_field_eq2(rho_n, rho_p)
+        protm = bsk.neutron_pairing_field(rho_n)
+        self.assertAlmostEqual(float(ref), float(protm), places=6)
+        
+    def test_proton_ref_pairing_field_eq2(self):
+        # Symmetric limit
         rho_n, rho_p = 0.05, 0.05
-        ref = bsk.ref_pairing_field_eq3(rho_n, rho_p)
+        ref = bsk.proton_ref_pairing_field_eq2(rho_n, rho_p)
         sym = bsk.symmetric_pairing_field(rho_n, rho_p)
         self.assertAlmostEqual(float(ref), float(sym), places=6)
+        # Pure neutron matter limit
+        rho_n, rho_p = 0.05, 0.00
+        ref   = bsk.proton_ref_pairing_field_eq2(rho_n, rho_p)
+        neum = bsk.neutron_pairing_field(rho_p) # Note: this is correct, the symmetry is assumed in this model
+        self.assertAlmostEqual(float(ref), float(neum), places=6)
+        # Pure proton matter limit
+        rho_n, rho_p = 0.00, 0.05
+        ref   = bsk.proton_ref_pairing_field_eq2(rho_n, rho_p)
+        protm = bsk.neutron_pairing_field(rho_p) # Note: this is correct, the symmetry is assumed in this model
+        self.assertAlmostEqual(float(ref), float(protm), places=6)
+
 
-    @unittest.skip("TODO(Adarsh): implement ref_pairing_field_eq6 (Eq. 6)")
-    def test_ref_pairing_field_eq6(self):
+    def test_neutron_ref_pairing_field_eq3(self):
+        # Symmetric limit
+        rho_n, rho_p = 0.05, 0.05
+        ref = bsk.neutron_ref_pairing_field_eq3(rho_n, rho_p)
+        sym = bsk.symmetric_pairing_field(rho_n, rho_p)
+        self.assertAlmostEqual(float(ref), float(sym), places=6)
+        # Pure neutron matter limit
+        rho_n, rho_p = 0.05, 0.00
+        ref = bsk.neutron_ref_pairing_field_eq3(rho_n, rho_p)
+        neum = bsk.neutron_pairing_field(rho_n)
+        self.assertAlmostEqual(float(ref), float(neum), places=6)
+        # Pure proton matter limit
+        rho_n, rho_p = 0.00, 0.05
+        ref = bsk.neutron_ref_pairing_field_eq3(rho_n, rho_p)
+        protm = bsk.neutron_pairing_field(rho_n) # Note: this is correct, the symmetry is assumed in this model
+        self.assertAlmostEqual(float(ref), float(protm), places=6)
+        
+    def test_proton_ref_pairing_field_eq3(self):
+        # Symmetric limit
+        rho_n, rho_p = 0.05, 0.05
+        ref = bsk.proton_ref_pairing_field_eq3(rho_n, rho_p)
+        sym = bsk.symmetric_pairing_field(rho_n, rho_p)
+        self.assertAlmostEqual(float(ref), float(sym), places=6)
+        # Pure neutron matter limit
+        rho_n, rho_p = 0.05, 0.00
+        ref = bsk.proton_ref_pairing_field_eq3(rho_n, rho_p)
+        neum = bsk.neutron_pairing_field(rho_p) # Note: this is correct, the symmetry is assumed in this model
+        self.assertAlmostEqual(float(ref), float(neum), places=6)
+        # Pure proton matter limit
+        rho_n, rho_p = 0.00, 0.05
+        ref = bsk.proton_ref_pairing_field_eq3(rho_n, rho_p)
+        protm = bsk.neutron_pairing_field(rho_p) # Note: this is correct, the symmetry is assumed in this model
+        self.assertAlmostEqual(float(ref), float(protm), places=6)
+
+    def test_neutron_ref_pairing_field_eq6(self):
+        # Symmetric limit
         rho_n, rho_p = 0.05, 0.05
-        ref = bsk.ref_pairing_field_eq6(rho_n, rho_p)
+        ref = bsk.neutron_ref_pairing_field_eq6(rho_n, rho_p)
         sym = bsk.symmetric_pairing_field(rho_n, rho_p)
         self.assertAlmostEqual(float(ref), float(sym), places=6)
+        # Pure neutron matter limit
+        rho_n, rho_p = 0.05, 0.00
+        ref = bsk.neutron_ref_pairing_field_eq6(rho_n, rho_p)
+        neum = bsk.neutron_pairing_field(rho_n)
+        self.assertAlmostEqual(float(ref), float(neum), places=6)
+        # Pure proton matter limit
+        rho_n, rho_p = 0.00, 0.05
+        ref = bsk.neutron_ref_pairing_field_eq6(rho_n, rho_p)
+        protm = bsk.neutron_pairing_field(rho_n) # Note: this is correct, the symmetry is assumed in this model
+        self.assertAlmostEqual(float(ref), float(protm), places=6)
+        
+    def test_proton_ref_pairing_field_eq6(self):
+        # Symmetric limit
+        rho_n, rho_p = 0.05, 0.05
+        ref = bsk.proton_ref_pairing_field_eq6(rho_n, rho_p)
+        sym = bsk.symmetric_pairing_field(rho_n, rho_p)
+        self.assertAlmostEqual(float(ref), float(sym), places=6)
+        # Pure neutron matter limit
+        rho_n, rho_p = 0.05, 0.00
+        ref = bsk.proton_ref_pairing_field_eq6(rho_n, rho_p)
+        neum = bsk.neutron_pairing_field(rho_p) # Note: this is correct, the symmetry is assumed in this model
+        self.assertAlmostEqual(float(ref), float(neum), places=6)
+        # Pure proton matter limit
+        rho_n, rho_p = 0.00, 0.05
+        ref = bsk.proton_ref_pairing_field_eq6(rho_n, rho_p)
+        protm = bsk.neutron_pairing_field(rho_p) # Note: this is correct, the symmetry is assumed in this model
+        self.assertAlmostEqual(float(ref), float(protm), places=6)
+        
+    def test_reference_fields_finite_scalar_and_array(self):
+        """All reference fields return finite values for scalar and array input."""
+        for func in (bsk.neutron_ref_pairing_field_eq2, bsk.neutron_ref_pairing_field_eq3,bsk.neutron_ref_pairing_field_eq6,
+                     bsk.proton_ref_pairing_field_eq2,bsk.proton_ref_pairing_field_eq3,bsk.proton_ref_pairing_field_eq6):
+            with self.subTest(func=func.__name__):
+                self.assertTrue(np.isfinite(func(0.05, 0.03)))
+                rho_n = np.array([0.02, 0.05, 0.08])
+                rho_p = np.array([0.01, 0.03, 0.04])
+                out = func(rho_n, rho_p)
+                self.assertEqual(out.shape, rho_n.shape)
+                self.assertTrue(np.all(np.isfinite(out)))
 
 
 if __name__ == '__main__':