diff --git a/mbtrack2/instability/__init__.py b/mbtrack2/instability/__init__.py
index f7a1516f07a5ba85219041527106425accf715ce..15fd68680eb6632db29e1eaf0a90369ddc7f4db3 100644
--- a/mbtrack2/instability/__init__.py
+++ b/mbtrack2/instability/__init__.py
@@ -7,6 +7,7 @@ from mbtrack2.instability.instabilities import (
     mbi_threshold,
     rwmbi_growth_rate,
     rwmbi_threshold,
+    transverse_gaussian_space_charge_tune_shift,
 )
 from mbtrack2.instability.ions import (
     fast_beam_ion,
diff --git a/mbtrack2/instability/instabilities.py b/mbtrack2/instability/instabilities.py
index 0a198eb98c90e9f55f16dc5a0805faa948889837..c67bc2205ce9b1de9642ed7c18f4210d128be374 100644
--- a/mbtrack2/instability/instabilities.py
+++ b/mbtrack2/instability/instabilities.py
@@ -454,3 +454,82 @@ def rwmbi_threshold(ring, beff, rho_material, plane='x'):
         (1-frac_tune) * omega0) / (2*c*Z0*rho_material))**0.5
 
     return Ith
+
+
+def transverse_gaussian_space_charge_tune_shift(ring, bunch_current, **kwargs):
+    """
+    Return the (maximum) transverse space charge tune shift for a Gaussian 
+    bunch in the linear approximation, see Eq.(1) of [1].
+
+    Parameters
+    ----------
+    ring : Synchrotron object
+        Ring parameters.
+    bunch_current : float
+        Bunch current in [A].
+    sigma_s : float, optional
+        RMS bunch length in [s].
+        Default is ring.sigma_0.
+    emit_x : float, optional
+        Horizontal emittance in [m.rad].
+        Default is ring.emit[0].
+    emit_y : float, optional
+        Vertical emittance in [m.rad].
+        Default is ring.emit[1].
+    use_lattice : bool, optional
+        If True, use beta fonctions along the lattice.
+        If False, local values of beta fonctions are used.
+        Default is ring.optics.use_local_values.
+    n_points : int, optional
+        Number of points in the lattice to be considered if use_lattice ==
+        True. Default is 1000.
+    sigma_delta : float, optional
+        Relative energy spread.
+        Default is ring.sigma_delta.
+    gamma : float, optional
+        Relativistic Lorentz gamma.
+        Default is ring.gamma.
+
+    Returns
+    -------
+    tune_shift : array of shape (2,)
+        Horizontal and vertical space charge tune shift.
+        
+    Reference
+    ---------
+    [1] : Antipov, S. A., Gubaidulin, V., Agapov, I., Garcia, E. C., & 
+    Gamelin, A. (2024). Space Charge and Future Light Sources. 
+    arXiv preprint arXiv:2409.08637.
+
+    """
+    sigma_s = kwargs.get('sigma_s', ring.sigma_0)
+    emit_x = kwargs.get('emit_x', ring.emit[0])
+    emit_y = kwargs.get('emit_y', ring.emit[1])
+    use_lattice = kwargs.get('use_lattice', not ring.optics.use_local_values)
+    sigma_delta = kwargs.get('sigma_delta', ring.sigma_delta)
+    gamma = kwargs.get('gamma', ring.gamma)
+    n_points = kwargs.get('n_points', 1000)
+    q = ring.particle.charge
+    m = ring.particle.mass
+    r_0 = 1 / (4*pi*epsilon_0) * q**2 / (m * c**2)
+    N = np.abs(bunch_current / ring.f0 / q)
+    sigma_z = sigma_s * c
+
+    if use_lattice:
+        s = np.linspace(0, ring.L, n_points)
+        beta = ring.optics.beta(s)
+        sig_x = (emit_x * beta[0] +
+                 ring.optics.dispersion(s)[0]**2 * sigma_delta**2)**0.5
+        sig_y = (emit_y * beta[1] +
+                 ring.optics.dispersion(s)[2]**2 * sigma_delta**2)**0.5
+        sig_xy = np.array([sig_x, sig_y])
+        return -r_0 * N / ((2 * pi)**1.5 * gamma**3 * sigma_z) * np.trapz(
+            beta / (sig_xy * (sig_y+sig_x)), s)
+    else:
+        beta = ring.optics.local_beta
+        sig_x = np.sqrt(beta[0] * emit_x)
+        sig_y = np.sqrt(beta[1] * emit_y)
+        sig_xy = np.array([sig_x, sig_y])
+        return -r_0 * N * ring.L / (
+            (2 * pi)**1.5 * gamma**3 * sigma_z) * beta / (sig_xy *
+                                                          (sig_y+sig_x))