diff --git a/stan/math/fwd.hpp b/stan/math/fwd.hpp
index 7f30955b0f2..3bab5be267f 100644
--- a/stan/math/fwd.hpp
+++ b/stan/math/fwd.hpp
@@ -11,6 +11,7 @@
 #include <stan/math/fwd/meta.hpp>
 #include <stan/math/fwd/fun.hpp>
 #include <stan/math/fwd/functor.hpp>
+#include <stan/math/fwd/prob.hpp>
 
 #include <stan/math/prim.hpp>
 
diff --git a/stan/math/fwd/prob.hpp b/stan/math/fwd/prob.hpp
new file mode 100644
index 00000000000..725acd8bb23
--- /dev/null
+++ b/stan/math/fwd/prob.hpp
@@ -0,0 +1,6 @@
+#ifndef STAN_MATH_FWD_PROB_HPP
+#define STAN_MATH_FWD_PROB_HPP
+
+#include <stan/math/fwd/prob/exponential_qf.hpp>
+
+#endif
diff --git a/stan/math/fwd/prob/exponential_qf.hpp b/stan/math/fwd/prob/exponential_qf.hpp
new file mode 100644
index 00000000000..c7f305d27f4
--- /dev/null
+++ b/stan/math/fwd/prob/exponential_qf.hpp
@@ -0,0 +1,67 @@
+#ifndef STAN_MATH_FWD_PROB_EXPONENTIAL_QF_HPP
+#define STAN_MATH_FWD_PROB_EXPONENTIAL_QF_HPP
+
+#include <stan/math/fwd/meta.hpp>
+#include <stan/math/fwd/core.hpp>
+#include <stan/math/prim/prob/exponential_qf.hpp>
+
+namespace stan {
+namespace math {
+
+template <typename Tp, typename Tbeta,
+          require_any_st_fvar<Tp, Tbeta>* = nullptr,
+          require_all_stan_scalar_t<Tp, Tbeta>* = nullptr>
+inline auto exponential_qf(const Tp& p, const Tbeta& beta) {
+  using fvar_t = return_type_t<Tp, Tbeta>;
+
+  auto p_val = value_of(p);
+  auto beta_val = value_of(beta);
+  fvar_t ret(exponential_qf(p_val, beta_val));
+
+  if (!is_constant<Tp>::value) {
+    ret.d_ += forward_as<fvar_t>(p).d_ * inv(beta_val - beta_val * p_val);
+  }
+  if (!is_constant<Tbeta>::value) {
+    ret.d_ += forward_as<fvar_t>(beta).d_ * log1m(p_val) / square(beta_val);
+  }
+
+  return ret;
+}
+
+template <typename Tp, typename Tbeta,
+          require_any_st_fvar<Tp, Tbeta>* = nullptr,
+          require_any_eigen_vector_t<Tp, Tbeta>* = nullptr>
+inline auto exponential_qf(const Tp& p, const Tbeta& beta) {
+  using vector_t
+      = plain_type_t<decltype(exponential_qf(value_of(p), value_of(beta)))>;
+  using fvar_t = return_type_t<Tp, Tbeta>;
+
+  auto p_val = value_of(p);
+  auto beta_val = value_of(beta);
+  vector_t ret = exponential_qf(p_val, beta_val);
+
+  auto p_array = as_array_or_scalar(p_val);
+  auto beta_array = as_array_or_scalar(beta_val);
+  vector_t d_ = vector_t::Zero(ret.rows(), ret.cols());
+
+  if (!is_constant<Tp>::value) {
+    as_array_or_scalar(d_)
+        += as_array_or_scalar(forward_as<promote_scalar_t<fvar_t, Tp>>(p).d())
+           * inv(beta_array - beta_array * p_array);
+  }
+  if (!is_constant<Tbeta>::value) {
+    as_array_or_scalar(d_)
+        += as_array_or_scalar(
+               forward_as<promote_scalar_t<fvar_t, Tbeta>>(beta).d())
+           * log1m(p_array) / square(beta_array);
+  }
+
+  return ret
+      .binaryExpr(d_, [&](const auto& val,
+                          const auto& deriv) { return fvar_t(val, deriv); })
+      .eval();
+}
+
+}  // namespace math
+}  // namespace stan
+#endif
diff --git a/stan/math/mix.hpp b/stan/math/mix.hpp
index 38b6a5e9c0e..63651abd76c 100644
--- a/stan/math/mix.hpp
+++ b/stan/math/mix.hpp
@@ -13,11 +13,13 @@
 #include <stan/math/fwd/meta.hpp>
 #include <stan/math/fwd/fun.hpp>
 #include <stan/math/fwd/functor.hpp>
+#include <stan/math/fwd/prob.hpp>
 
 #include <stan/math/rev/core.hpp>
 #include <stan/math/rev/meta.hpp>
 #include <stan/math/rev/fun.hpp>
 #include <stan/math/rev/functor.hpp>
+#include <stan/math/rev/prob.hpp>
 
 #include <stan/math/prim.hpp>
 
diff --git a/stan/math/prim/prob.hpp b/stan/math/prim/prob.hpp
index 33584411ec3..77cf4ed8ee0 100644
--- a/stan/math/prim/prob.hpp
+++ b/stan/math/prim/prob.hpp
@@ -106,6 +106,7 @@
 #include <stan/math/prim/prob/exponential_lcdf.hpp>
 #include <stan/math/prim/prob/exponential_log.hpp>
 #include <stan/math/prim/prob/exponential_lpdf.hpp>
+#include <stan/math/prim/prob/exponential_qf.hpp>
 #include <stan/math/prim/prob/exponential_rng.hpp>
 #include <stan/math/prim/prob/frechet_ccdf_log.hpp>
 #include <stan/math/prim/prob/frechet_cdf.hpp>
diff --git a/stan/math/prim/prob/exponential_qf.hpp b/stan/math/prim/prob/exponential_qf.hpp
new file mode 100644
index 00000000000..7d4fa96e0fd
--- /dev/null
+++ b/stan/math/prim/prob/exponential_qf.hpp
@@ -0,0 +1,70 @@
+#ifndef STAN_MATH_PRIM_PROB_EXPONENTIAL_QF_HPP
+#define STAN_MATH_PRIM_PROB_EXPONENTIAL_QF_HPP
+
+#include <stan/math/prim/meta.hpp>
+#include <stan/math/prim/err.hpp>
+#include <stan/math/prim/fun/Eigen.hpp>
+#include <vector>
+
+namespace stan {
+namespace math {
+
+/** \ingroup prob_dists
+ * The quantile of an exponential density for p with the specified
+ * inverse scale parameter.
+ * Inverse scale parameter must be greater than 0.
+ * p must be bounded by 0 and 1.
+ *
+ \f[
+   q = -\log(1 - p) / \beta
+ \f]
+ *
+ * @tparam Tp type of probability input
+ * @tparam Tbeta type of inverse scale
+ * @param p A scalar variable.
+ * @param beta Inverse scale parameter.
+ * @throw std::domain_error if beta is not greater than 0.
+ * @throw std::domain_error if y is not greater than or equal to 0.
+ */
+template <typename Tp, typename Tbeta,
+          require_all_arithmetic_t<Tp, Tbeta>* = nullptr>
+inline auto exponential_qf(const Tp& p, const Tbeta& beta) {
+  static const char* function = "exponential_qf";
+  check_positive_finite(function, "Inverse scale parameter", beta);
+  check_bounded(function, "Probability parameter", p, 0, 1);
+  return -log1p(-p) / beta;
+}
+
+/** \ingroup prob_dists
+ * The quantile of an exponential density for p with the specified
+ * inverse scale parameter.
+ * Inverse scale parameter must be greater than 0.
+ * p must be bounded by 0 and 1.
+ *
+ * Specialisation for use where any input is an Eigen vector
+ *
+ * @tparam Tp type of probability input
+ * @tparam Tbeta type of inverse scale
+ * @param p A vector or scalar of probabilities.
+ * @param beta Vector or scalar of inverse scale parameter.
+ * @throw std::domain_error if beta is not greater than 0.
+ * @throw std::domain_error if y is not greater than or equal to 0.
+ */
+template <typename Tp, typename Tbeta,
+          require_all_st_arithmetic<Tp, Tbeta>* = nullptr,
+          require_any_eigen_vector_t<Tp, Tbeta>* = nullptr>
+inline auto exponential_qf(const Tp& p, const Tbeta& beta) {
+  static const char* function = "exponential_qf";
+  ref_type_t<Tp> p_ref = p;
+  ref_type_t<Tbeta> beta_ref = beta;
+  check_positive_finite(function, "Inverse scale parameter", beta_ref);
+  check_bounded(function, "Probability parameter", p_ref, 0, 1);
+  return (-log1p(-as_array_or_scalar(p_ref)) / as_array_or_scalar(beta_ref))
+      .matrix()
+      .eval();
+}
+
+}  // namespace math
+}  // namespace stan
+
+#endif
diff --git a/stan/math/rev.hpp b/stan/math/rev.hpp
index 866deb448d3..20d5a477846 100644
--- a/stan/math/rev.hpp
+++ b/stan/math/rev.hpp
@@ -9,6 +9,7 @@
 #include <stan/math/rev/meta.hpp>
 #include <stan/math/rev/fun.hpp>
 #include <stan/math/rev/functor.hpp>
+#include <stan/math/rev/prob.hpp>
 
 #include <stan/math/prim.hpp>
 
diff --git a/stan/math/rev/prob.hpp b/stan/math/rev/prob.hpp
new file mode 100644
index 00000000000..48156af85ed
--- /dev/null
+++ b/stan/math/rev/prob.hpp
@@ -0,0 +1,6 @@
+#ifndef STAN_MATH_REV_PROB_HPP
+#define STAN_MATH_REV_PROB_HPP
+
+#include <stan/math/rev/prob/exponential_qf.hpp>
+
+#endif
diff --git a/stan/math/rev/prob/exponential_qf.hpp b/stan/math/rev/prob/exponential_qf.hpp
new file mode 100644
index 00000000000..ff2d591ee00
--- /dev/null
+++ b/stan/math/rev/prob/exponential_qf.hpp
@@ -0,0 +1,57 @@
+#ifndef STAN_MATH_REV_PROB_EXPONENTIAL_QF_HPP
+#define STAN_MATH_REV_PROB_EXPONENTIAL_QF_HPP
+
+#include <stan/math/rev/meta.hpp>
+#include <stan/math/rev/core.hpp>
+#include <stan/math/prim/prob/exponential_qf.hpp>
+
+namespace stan {
+namespace math {
+
+/** \ingroup prob_dists
+ * The quantile of an exponential density for p with the specified
+ * inverse scale parameter.
+ * Inverse scale parameter must be greater than 0.
+ * p must be bounded by 0 and 1.
+ *
+ \f[
+   \frac{\partial }{\partial p} = \frac{1}{\beta - \beta p} \\
+   \frac{\partial }{\partial \beta} = \frac{\log(1 - p)}{\beta^2}
+ \f]
+ *
+ * @tparam Tp type of probability input
+ * @tparam Tbeta type of inverse scale
+ * @param p A scalar variable.
+ * @param beta Inverse scale parameter.
+ * @throw std::domain_error if beta is not greater than 0.
+ * @throw std::domain_error if y is not greater than or equal to 0.
+ */
+template <typename Tp, typename Tbeta, require_any_st_var<Tp, Tbeta>* = nullptr>
+inline auto exponential_qf(const Tp& p, const Tbeta& beta) {
+  using inner_ret_type = decltype(exponential_qf(value_of(p), value_of(beta)));
+  using ret_type = return_var_matrix_t<inner_ret_type, Tp, Tbeta>;
+
+  arena_t<Tp> arena_p = p;
+  arena_t<Tbeta> arena_beta = beta;
+  arena_t<ret_type> ret
+      = exponential_qf(value_of(arena_p), value_of(arena_beta));
+  reverse_pass_callback([ret, arena_p, arena_beta]() mutable {
+    decltype(auto) p_val = as_array_or_scalar(value_of(arena_p));
+    decltype(auto) beta_val = as_array_or_scalar(value_of(arena_beta));
+    decltype(auto) ret_adj = as_array_or_scalar(ret.adj());
+    if (!is_constant<Tp>::value) {
+      as_array_or_scalar(forward_as<promote_scalar_t<var, Tp>>(arena_p).adj())
+          += ret_adj * inv(beta_val - beta_val * p_val);
+    }
+    if (!is_constant<Tbeta>::value) {
+      as_array_or_scalar(
+          forward_as<promote_scalar_t<var, Tbeta>>(arena_beta).adj())
+          += ret_adj * log1m(p_val) / square(beta_val);
+    }
+  });
+  return ret_type(ret);
+}
+
+}  // namespace math
+}  // namespace stan
+#endif
diff --git a/test/unit/math/mix/prob/exponential_qf_test.cpp b/test/unit/math/mix/prob/exponential_qf_test.cpp
new file mode 100644
index 00000000000..8fe52f4749e
--- /dev/null
+++ b/test/unit/math/mix/prob/exponential_qf_test.cpp
@@ -0,0 +1,17 @@
+#include <stan/math/mix.hpp>
+#include <test/unit/math/test_ad.hpp>
+
+TEST(mathMixProb, exponential_qf) {
+  auto f = [](const auto& p, const auto& beta) {
+    return stan::math::exponential_qf(p, beta);
+  };
+
+  Eigen::VectorXd p(2);
+  p << 0.2, 0.6;
+
+  Eigen::VectorXd beta(2);
+  beta << 16, 2;
+
+  stan::test::expect_ad(f, p, beta);
+  stan::test::expect_ad(f, 0.1, 8);
+}
diff --git a/test/unit/math/prim/prob/exponential_qf_test.cpp b/test/unit/math/prim/prob/exponential_qf_test.cpp
new file mode 100644
index 00000000000..2ceca7d76c0
--- /dev/null
+++ b/test/unit/math/prim/prob/exponential_qf_test.cpp
@@ -0,0 +1,53 @@
+#include <stan/math/prim.hpp>
+#include <gtest/gtest.h>
+#include <test/unit/util.hpp>
+
+TEST(ProbExponential, quantile_values) {
+  using stan::math::exponential_qf;
+
+  EXPECT_FLOAT_EQ(exponential_qf(0.1, 8), 0.0131700644572);
+  EXPECT_FLOAT_EQ(exponential_qf(0.4, 7), 0.0729750891094);
+  EXPECT_FLOAT_EQ(exponential_qf(0.1, 7), 0.0150515022368);
+
+  Eigen::VectorXd p(2);
+  p << 0.2, 0.6;
+
+  Eigen::VectorXd beta(2);
+  beta << 16, 2;
+
+  Eigen::VectorXd res(2);
+  res << 0.0318776501877, 0.1308986759820;
+
+  EXPECT_MATRIX_FLOAT_EQ(exponential_qf(p, 7), res);
+
+  res << 0.0222921839962, 0.1783374719694;
+  EXPECT_MATRIX_FLOAT_EQ(exponential_qf(0.3, beta), res);
+
+  res << 0.0139464719571, 0.4581453659371;
+  EXPECT_MATRIX_FLOAT_EQ(exponential_qf(p, beta), res);
+}
+
+TEST(ProbExponential, quantile_throws) {
+  using stan::math::exponential_qf;
+
+  EXPECT_THROW(exponential_qf(2, 8), std::domain_error);
+  EXPECT_THROW(exponential_qf(0.4, -7), std::domain_error);
+
+  Eigen::VectorXd p(2);
+  p << 0.2, 0.6;
+
+  Eigen::VectorXd p_invalid(2);
+  p_invalid << 0.2, 2.2;
+
+  Eigen::VectorXd beta(2);
+  beta << 16, 2;
+
+  Eigen::VectorXd beta_invalid(2);
+  beta_invalid << 16, -8;
+
+  EXPECT_THROW(exponential_qf(p_invalid, 8), std::domain_error);
+  EXPECT_THROW(exponential_qf(p_invalid, beta), std::domain_error);
+  EXPECT_THROW(exponential_qf(0.4, beta_invalid), std::domain_error);
+  EXPECT_THROW(exponential_qf(p, beta_invalid), std::domain_error);
+  EXPECT_THROW(exponential_qf(p_invalid, beta_invalid), std::domain_error);
+}