diff --git a/Grid/qcd/modules/ObservableModules.h b/Grid/qcd/modules/ObservableModules.h
index 87fcbb92..cda86185 100644
--- a/Grid/qcd/modules/ObservableModules.h
+++ b/Grid/qcd/modules/ObservableModules.h
@@ -103,6 +103,18 @@ class PolyakovMod: public ObservableModule<PolyakovLogger<Impl>, NoParameters>{
   PolyakovMod(): ObsBase(NoParameters()){}
 };
 
+template < class Impl >
+class SpatialPolyakovMod: public ObservableModule<SpatialPolyakovLogger<Impl>, NoParameters>{
+  typedef ObservableModule<SpatialPolyakovLogger<Impl>, NoParameters> ObsBase;
+  using ObsBase::ObsBase; // for constructors
+
+  // acquire resource
+  virtual void initialize(){
+    this->ObservablePtr.reset(new SpatialPolyakovLogger<Impl>());
+  }
+  public:
+  SpatialPolyakovMod(): ObsBase(NoParameters()){}
+};
 
 template < class Impl >
 class TopologicalChargeMod: public ObservableModule<TopologicalCharge<Impl>, TopologyObsParameters>{
diff --git a/Grid/qcd/observables/polyakov_loop.h b/Grid/qcd/observables/polyakov_loop.h
index 0b59f549..57812ff6 100644
--- a/Grid/qcd/observables/polyakov_loop.h
+++ b/Grid/qcd/observables/polyakov_loop.h
@@ -2,11 +2,12 @@
 
 Grid physics library, www.github.com/paboyle/Grid
 
-Source file: ./lib/qcd/modules/polyakov_line.h
+Source file: ./Grid/qcd/observables/polyakov_loop.h
 
-Copyright (C) 2017
+Copyright (C) 2025
 
 Author: David Preti <david.preti@csic.es>
+Author: Alexis Verney-Provatas <2414441@swansea.ac.uk>
 
 This program is free software; you can redistribute it and/or modify
 it under the terms of the GNU General Public License as published by
@@ -60,4 +61,43 @@ class PolyakovLogger : public HmcObservable<typename Impl::Field> {
   }
 };
 
+template <class Impl>
+class SpatialPolyakovLogger : public HmcObservable<typename Impl::Field> {
+ public:
+    // here forces the Impl to be of gauge fields
+    // if not the compiler will complain
+    INHERIT_GIMPL_TYPES(Impl);
+
+     // necessary for HmcObservable compatibility
+    typedef typename Impl::Field Field;
+
+    void TrajectoryComplete(int traj,
+                            Field &U,
+                            GridSerialRNG &sRNG,
+                            GridParallelRNG &pRNG) {
+
+    // Save current numerical output precision
+    int def_prec = std::cout.precision();
+
+    // Assume that the dimensions are D=3+1
+    int Ndim = 3;
+    ComplexD polyakov;
+   
+    // Iterate over the spatial directions and print the average spatial polyakov loop
+    // over them 
+    for (int idx=0; idx<Ndim; idx++) {
+        polyakov = WilsonLoops<Impl>::avgPolyakovLoop(U, idx);
+    
+        std::cout << GridLogMessage
+            << std::setprecision(std::numeric_limits<Real>::digits10 + 1)
+            << "Polyakov Loop in the " << idx << " spatial direction : [ " << traj << " ] "<< polyakov << std::endl;
+
+    }
+
+    // Return to original output precision
+    std::cout.precision(def_prec);
+
+  }
+};
+
 NAMESPACE_END(Grid);
diff --git a/Grid/qcd/utils/WilsonLoops.h b/Grid/qcd/utils/WilsonLoops.h
index 23ad237c..8091cbc8 100644
--- a/Grid/qcd/utils/WilsonLoops.h
+++ b/Grid/qcd/utils/WilsonLoops.h
@@ -177,25 +177,43 @@ public:
   }
 
   //////////////////////////////////////////////////
-  // average over all x,y,z the temporal loop
+  // average Polyakov loop in mu direction over all directions != mu
   //////////////////////////////////////////////////
-  static ComplexD avgPolyakovLoop(const GaugeField &Umu) {  //assume Nd=4
-    GaugeMat Ut(Umu.Grid()), P(Umu.Grid());
+  static ComplexD avgPolyakovLoop(const GaugeField &Umu, const int mu) {  //assume Nd=4
+    
+    // Protect against bad value of mu [0, 3]
+    if ((mu < 0 ) || (mu > 3)) {
+      std::cout << GridLogError << "Index is not an integer inclusively between 0 and 3." << std::endl;
+      exit(1);
+    }
+
+    // U_loop is U_{mu}
+    GaugeMat U_loop(Umu.Grid()), P(Umu.Grid());
     ComplexD out;
     int T = Umu.Grid()->GlobalDimensions()[3];
     int X = Umu.Grid()->GlobalDimensions()[0];
     int Y = Umu.Grid()->GlobalDimensions()[1];
     int Z = Umu.Grid()->GlobalDimensions()[2];
 
-    Ut = peekLorentz(Umu,3); //Select temporal direction
-    P = Ut;
-    for (int t=1;t<T;t++){ 
-      P = Gimpl::CovShiftForward(Ut,3,P);
+    // Number of sites in mu direction
+    int N_mu = Umu.Grid()->GlobalDimensions()[mu];
+
+    U_loop = peekLorentz(Umu, mu); //Select direction
+    P = U_loop;
+    for (int t=1;t<N_mu;t++){ 
+      P = Gimpl::CovShiftForward(U_loop,mu,P);
     }
    RealD norm = 1.0/(Nc*X*Y*Z*T);
    out = sum(trace(P))*norm;
    return out;   
-}
+  }  
+
+  /////////////////////////////////////////////////
+  // overload for temporal Polyakov loop
+  /////////////////////////////////////////////////
+  static ComplexD avgPolyakovLoop(const GaugeField &Umu) { 
+    return avgPolyakovLoop(Umu, 3);
+  }
 
   //////////////////////////////////////////////////
   // average over traced single links