diff --git a/src/language/modules/LocalFSIModule.cpp b/src/language/modules/LocalFSIModule.cpp
index 878a93faae496b6c16a74f05a633b066c5eb6c77..4c46e950fd7084a7672e2d3c42d270d2e0bbbb93 100644
--- a/src/language/modules/LocalFSIModule.cpp
+++ b/src/language/modules/LocalFSIModule.cpp
@@ -1077,6 +1077,22 @@ LocalFSIModule::LocalFSIModule()
                               }
 
                               ));
+                            
+  this->_addBuiltinFunction("acoustic_dt_order3", std::function(
+
+                                             [](const std::shared_ptr<const DiscreteFunctionVariant>& c) -> double {
+                                               return acoustic_dt_order3(c);
+                                             }
+
+                                             ));
+
+  this->_addBuiltinFunction("acoustic_dt_order4", std::function(
+
+                                             [](const std::shared_ptr<const DiscreteFunctionVariant>& c) -> double {
+                                               return acoustic_dt_order4(c);
+                                             }
+
+                                             ));
 }
 
 void
diff --git a/src/scheme/HighOrderAcousticSolver.cpp b/src/scheme/HighOrderAcousticSolver.cpp
index 8db1716542ed119aba814acd5073d00729884619..dc7c70a5c2daa5f60ff501188b8bc72e57764833 100644
--- a/src/scheme/HighOrderAcousticSolver.cpp
+++ b/src/scheme/HighOrderAcousticSolver.cpp
@@ -30,6 +30,58 @@
 #include <variant>
 #include <vector>
 
+double 
+acoustic_dt_order3(const std::shared_ptr<const DiscreteFunctionVariant>& c_v)
+{
+  const auto& c = c_v->get<DiscreteFunctionP0<const double>>();
+
+  return std::visit(
+    [&](auto&& p_mesh) -> double {
+      const auto& mesh = *p_mesh;
+
+      using MeshType = decltype(mesh);
+      if constexpr (is_polygonal_mesh_v<MeshType>) {
+        const auto Vj = MeshDataManager::instance().getMeshData(mesh).Vj();
+        const auto Sj = MeshDataManager::instance().getMeshData(mesh).sumOverRLjr();
+
+        CellValue<double> local_dt{mesh.connectivity()};
+        parallel_for(
+          mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { local_dt[j] = 2 * pow(Vj[j],3/2) / (Sj[j] * c[j]); });
+
+        return min(local_dt);
+      } else {
+        throw NormalError("unexpected mesh type");
+      }
+    },
+    c.meshVariant()->variant());
+}
+
+double 
+acoustic_dt_order4(const std::shared_ptr<const DiscreteFunctionVariant>& c_v)
+{
+  const auto& c = c_v->get<DiscreteFunctionP0<const double>>();
+
+  return std::visit(
+    [&](auto&& p_mesh) -> double {
+      const auto& mesh = *p_mesh;
+
+      using MeshType = decltype(mesh);
+      if constexpr (is_polygonal_mesh_v<MeshType>) {
+        const auto Vj = MeshDataManager::instance().getMeshData(mesh).Vj();
+        const auto Sj = MeshDataManager::instance().getMeshData(mesh).sumOverRLjr();
+
+        CellValue<double> local_dt{mesh.connectivity()};
+        parallel_for(
+          mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { local_dt[j] = 2 * Vj[j] * Vj[j] / (Sj[j] * c[j]); });
+
+        return min(local_dt);
+      } else {
+        throw NormalError("unexpected mesh type");
+      }
+    },
+    c.meshVariant()->variant());
+}
+
 template <MeshConcept MeshType>
 class HighOrderAcousticSolverHandler::HighOrderAcousticSolver final : public HighOrderAcousticSolverHandler::IHighOrderAcousticSolver
 {
diff --git a/src/scheme/HighOrderAcousticSolver.hpp b/src/scheme/HighOrderAcousticSolver.hpp
index 8545750ccf4ced9e197427387df42c6222afefd0..a85522cb19c763cf147775245c1775fd345edd32 100644
--- a/src/scheme/HighOrderAcousticSolver.hpp
+++ b/src/scheme/HighOrderAcousticSolver.hpp
@@ -13,6 +13,10 @@ class MeshVariant;
 class ItemValueVariant;
 class SubItemValuePerItemVariant;
 
+double acoustic_dt_order3(const std::shared_ptr<const DiscreteFunctionVariant>& c_v);
+
+double acoustic_dt_order4(const std::shared_ptr<const DiscreteFunctionVariant>& c_v);
+
 class HighOrderAcousticSolverHandler
 {
  public: