diff --git a/stan/math/prim/constraint.hpp b/stan/math/prim/constraint.hpp
index 169fd8f2e31..5c494c83214 100644
--- a/stan/math/prim/constraint.hpp
+++ b/stan/math/prim/constraint.hpp
@@ -37,6 +37,8 @@
 #include <stan/math/prim/constraint/stochastic_column_free.hpp>
 #include <stan/math/prim/constraint/stochastic_row_constrain.hpp>
 #include <stan/math/prim/constraint/stochastic_row_free.hpp>
+#include <stan/math/prim/constraint/sum_to_zero_constrain.hpp>
+#include <stan/math/prim/constraint/sum_to_zero_free.hpp>
 #include <stan/math/prim/constraint/ub_constrain.hpp>
 #include <stan/math/prim/constraint/ub_free.hpp>
 #include <stan/math/prim/constraint/unit_vector_constrain.hpp>
diff --git a/stan/math/prim/constraint/sum_to_zero_constrain.hpp b/stan/math/prim/constraint/sum_to_zero_constrain.hpp
new file mode 100644
index 00000000000..c3979720c11
--- /dev/null
+++ b/stan/math/prim/constraint/sum_to_zero_constrain.hpp
@@ -0,0 +1,166 @@
+#ifndef STAN_MATH_PRIM_CONSTRAINT_SUM_TO_ZERO_CONSTRAIN_HPP
+#define STAN_MATH_PRIM_CONSTRAINT_SUM_TO_ZERO_CONSTRAIN_HPP
+
+#include <stan/math/prim/meta.hpp>
+#include <stan/math/prim/fun/Eigen.hpp>
+#include <stan/math/prim/fun/to_ref.hpp>
+#include <stan/math/prim/fun/inv_sqrt.hpp>
+#include <stan/math/prim/functor/apply_vector_unary.hpp>
+#include <cmath>
+
+namespace stan {
+namespace math {
+
+/**
+ * Return a vector with sum zero corresponding to the specified
+ * free vector.
+ *
+ * The sum-to-zero transform is defined using a modified version of the
+ * the inverse of the isometric log ratio transform (ILR).
+ * See:
+ * Egozcue, Juan Jose; Pawlowsky-Glahn, Vera; Mateu-Figueras, Gloria;
+ * Barcelo-Vidal, Carles (2003), "Isometric logratio transformations for
+ * compositional data analysis", Mathematical Geology, 35 (3): 279–300,
+ * doi:10.1023/A:1023818214614, S2CID 122844634
+ *
+ * This implementation is closer to the description of the same using "pivot
+ * coordinates" in
+ * Filzmoser, P., Hron, K., Templ, M. (2018). Geometrical Properties of
+ * Compositional Data. In: Applied Compositional Data Analysis. Springer Series
+ * in Statistics. Springer, Cham. https://doi.org/10.1007/978-3-319-96422-5_3
+ *
+ * This is a linear transform, with no Jacobian.
+ *
+ * @tparam Vec type of the vector
+ * @param y Free vector input of dimensionality K - 1.
+ * @return Zero-sum vector of dimensionality K.
+ */
+template <typename Vec, require_eigen_col_vector_t<Vec>* = nullptr,
+          require_not_st_var<Vec>* = nullptr>
+inline plain_type_t<Vec> sum_to_zero_constrain(const Vec& y) {
+  const auto N = y.size();
+
+  plain_type_t<Vec> z = Eigen::VectorXd::Zero(N + 1);
+  if (unlikely(N == 0)) {
+    return z;
+  }
+
+  auto&& y_ref = to_ref(y);
+
+  value_type_t<Vec> sum_w(0);
+  for (int i = N; i > 0; --i) {
+    double n = static_cast<double>(i);
+    auto w = y_ref(i - 1) * inv_sqrt(n * (n + 1));
+    sum_w += w;
+
+    z.coeffRef(i - 1) += sum_w;
+    z.coeffRef(i) -= w * n;
+  }
+
+  return z;
+}
+
+/**
+ * Return a vector with sum zero corresponding to the specified
+ * free vector.
+ *
+ * The sum-to-zero transform is defined using a modified version of the
+ * the inverse of the isometric log ratio transform (ILR).
+ * See:
+ * Egozcue, Juan Jose; Pawlowsky-Glahn, Vera; Mateu-Figueras, Gloria;
+ * Barcelo-Vidal, Carles (2003), "Isometric logratio transformations for
+ * compositional data analysis", Mathematical Geology, 35 (3): 279–300,
+ * doi:10.1023/A:1023818214614, S2CID 122844634
+ *
+ * This implementation is closer to the description of the same using "pivot
+ * coordinates" in
+ * Filzmoser, P., Hron, K., Templ, M. (2018). Geometrical Properties of
+ * Compositional Data. In: Applied Compositional Data Analysis. Springer Series
+ * in Statistics. Springer, Cham. https://doi.org/10.1007/978-3-319-96422-5_3
+ *
+ * This is a linear transform, with no Jacobian.
+ *
+ * @tparam Vec type of the vector
+ * @param y Free vector input of dimensionality K - 1.
+ * @param lp unused
+ * @return Zero-sum vector of dimensionality K.
+ */
+template <typename Vec, require_eigen_col_vector_t<Vec>* = nullptr,
+          require_not_st_var<Vec>* = nullptr>
+inline plain_type_t<Vec> sum_to_zero_constrain(const Vec& y,
+                                               value_type_t<Vec>& lp) {
+  return sum_to_zero_constrain(y);
+}
+
+/**
+ * Return a vector with sum zero corresponding to the specified
+ * free vector.
+ *
+ * The sum-to-zero transform is defined using a modified version of
+ * the inverse of the isometric log ratio transform (ILR).
+ * See:
+ * Egozcue, Juan Jose; Pawlowsky-Glahn, Vera; Mateu-Figueras, Gloria;
+ * Barcelo-Vidal, Carles (2003), "Isometric logratio transformations for
+ * compositional data analysis", Mathematical Geology, 35 (3): 279–300,
+ * doi:10.1023/A:1023818214614, S2CID 122844634
+ *
+ * This implementation is closer to the description of the same using "pivot
+ * coordinates" in
+ * Filzmoser, P., Hron, K., Templ, M. (2018). Geometrical Properties of
+ * Compositional Data. In: Applied Compositional Data Analysis. Springer Series
+ * in Statistics. Springer, Cham. https://doi.org/10.1007/978-3-319-96422-5_3
+ *
+ * This is a linear transform, with no Jacobian.
+ *
+ * @tparam Jacobian unused
+ * @tparam Vec A type inheriting from `Eigen::DenseBase` or a `var_value` with
+ *  inner type inheriting from `Eigen::DenseBase` with compile time dynamic rows
+ *  and 1 column
+ * @param[in] y free vector
+ * @param[in, out] lp unused
+ * @return Zero-sum vector of dimensionality one greater than `y`
+ */
+template <bool Jacobian, typename Vec, require_not_std_vector_t<Vec>* = nullptr>
+inline plain_type_t<Vec> sum_to_zero_constrain(const Vec& y,
+                                               return_type_t<Vec>& lp) {
+  return sum_to_zero_constrain(y);
+}
+
+/**
+ * Return a vector with sum zero corresponding to the specified
+ * free vector.
+ *
+ * The sum-to-zero transform is defined using a modified version of
+ * the inverse of the isometric log ratio transform (ILR).
+ * See:
+ * Egozcue, Juan Jose; Pawlowsky-Glahn, Vera; Mateu-Figueras, Gloria;
+ * Barcelo-Vidal, Carles (2003), "Isometric logratio transformations for
+ * compositional data analysis", Mathematical Geology, 35 (3): 279–300,
+ * doi:10.1023/A:1023818214614, S2CID 122844634
+ *
+ * This implementation is closer to the description of the same using "pivot
+ * coordinates" in
+ * Filzmoser, P., Hron, K., Templ, M. (2018). Geometrical Properties of
+ * Compositional Data. In: Applied Compositional Data Analysis. Springer Series
+ * in Statistics. Springer, Cham. https://doi.org/10.1007/978-3-319-96422-5_3
+ *
+ * This is a linear transform, with no Jacobian.
+ *
+ * @tparam Jacobian unused
+ * @tparam Vec A standard vector with inner type inheriting from
+ * `Eigen::DenseBase` or a `var_value` with inner type inheriting from
+ * `Eigen::DenseBase` with compile time dynamic rows and 1 column
+ * @param[in] y free vector
+ * @param[in, out] lp unused
+ * @return Zero-sum vectors of dimensionality one greater than `y`
+ */
+template <bool Jacobian, typename T, require_std_vector_t<T>* = nullptr>
+inline auto sum_to_zero_constrain(const T& y, return_type_t<T>& lp) {
+  return apply_vector_unary<T>::apply(
+      y, [](auto&& v) { return sum_to_zero_constrain(v); });
+}
+
+}  // namespace math
+}  // namespace stan
+
+#endif
diff --git a/stan/math/prim/constraint/sum_to_zero_free.hpp b/stan/math/prim/constraint/sum_to_zero_free.hpp
new file mode 100644
index 00000000000..650970f992d
--- /dev/null
+++ b/stan/math/prim/constraint/sum_to_zero_free.hpp
@@ -0,0 +1,80 @@
+#ifndef STAN_MATH_PRIM_CONSTRAINT_SUM_TO_ZERO_FREE_HPP
+#define STAN_MATH_PRIM_CONSTRAINT_SUM_TO_ZERO_FREE_HPP
+
+#include <stan/math/prim/meta.hpp>
+#include <stan/math/prim/err.hpp>
+#include <stan/math/prim/fun/Eigen.hpp>
+#include <stan/math/prim/fun/to_ref.hpp>
+#include <stan/math/prim/fun/sqrt.hpp>
+#include <stan/math/prim/functor/apply_vector_unary.hpp>
+#include <cmath>
+
+namespace stan {
+namespace math {
+
+/**
+ * Return an unconstrained vector.
+ *
+ * The sum-to-zero transform is defined using a modified version of the
+ * isometric log ratio transform (ILR).
+ * See:
+ * Egozcue, Juan Jose; Pawlowsky-Glahn, Vera; Mateu-Figueras, Gloria;
+ * Barcelo-Vidal, Carles (2003), "Isometric logratio transformations for
+ * compositional data analysis", Mathematical Geology, 35 (3): 279–300,
+ * doi:10.1023/A:1023818214614, S2CID 122844634
+ *
+ * This implementation is closer to the description of the same using "pivot
+ * coordinates" in
+ * Filzmoser, P., Hron, K., Templ, M. (2018). Geometrical Properties of
+ * Compositional Data. In: Applied Compositional Data Analysis. Springer Series
+ * in Statistics. Springer, Cham. https://doi.org/10.1007/978-3-319-96422-5_3
+ *
+ * @tparam ColVec a column vector type
+ * @param z Vector of length K.
+ * @return Free vector of length (K-1).
+ * @throw std::domain_error if z does not sum to zero
+ */
+template <typename Vec, require_eigen_vector_t<Vec>* = nullptr>
+inline plain_type_t<Vec> sum_to_zero_free(const Vec& z) {
+  const auto& z_ref = to_ref(z);
+  check_sum_to_zero("stan::math::sum_to_zero_free", "sum_to_zero variable",
+                    z_ref);
+
+  const auto N = z.size() - 1;
+
+  plain_type_t<Vec> y = Eigen::VectorXd::Zero(N);
+  if (unlikely(N == 0)) {
+    return y;
+  }
+
+  y.coeffRef(N - 1) = -z_ref(N) * sqrt(N * (N + 1)) / N;
+
+  value_type_t<Vec> sum_w(0);
+
+  for (int i = N - 2; i >= 0; --i) {
+    double n = static_cast<double>(i + 1);
+    auto w = y(i + 1) / sqrt((n + 1) * (n + 2));
+    sum_w += w;
+    y.coeffRef(i) = (sum_w - z_ref(i + 1)) * sqrt(n * (n + 1)) / n;
+  }
+
+  return y;
+}
+
+/**
+ * Overload of `sum_to_zero_free()` to untransform each Eigen vector
+ * in a standard vector.
+ * @tparam T A standard vector with with a `value_type` which inherits from
+ *  `Eigen::MatrixBase` with compile time rows or columns equal to 1.
+ * @param z The standard vector to untransform.
+ */
+template <typename T, require_std_vector_t<T>* = nullptr>
+auto sum_to_zero_free(const T& z) {
+  return apply_vector_unary<T>::apply(
+      z, [](auto&& v) { return sum_to_zero_free(v); });
+}
+
+}  // namespace math
+}  // namespace stan
+
+#endif
diff --git a/stan/math/prim/err.hpp b/stan/math/prim/err.hpp
index c78344c95d8..72d6c733e39 100644
--- a/stan/math/prim/err.hpp
+++ b/stan/math/prim/err.hpp
@@ -41,6 +41,7 @@
 #include <stan/math/prim/err/check_std_vector_index.hpp>
 #include <stan/math/prim/err/check_stochastic_column.hpp>
 #include <stan/math/prim/err/check_stochastic_row.hpp>
+#include <stan/math/prim/err/check_sum_to_zero.hpp>
 #include <stan/math/prim/err/check_symmetric.hpp>
 #include <stan/math/prim/err/check_unit_vector.hpp>
 #include <stan/math/prim/err/check_vector.hpp>
diff --git a/stan/math/prim/err/check_sum_to_zero.hpp b/stan/math/prim/err/check_sum_to_zero.hpp
new file mode 100644
index 00000000000..b384f5cc778
--- /dev/null
+++ b/stan/math/prim/err/check_sum_to_zero.hpp
@@ -0,0 +1,72 @@
+#ifndef STAN_MATH_PRIM_ERR_CHECK_SUM_TO_ZERO_HPP
+#define STAN_MATH_PRIM_ERR_CHECK_SUM_TO_ZERO_HPP
+
+#include <stan/math/prim/fun/Eigen.hpp>
+#include <stan/math/prim/meta.hpp>
+#include <stan/math/prim/err/constraint_tolerance.hpp>
+#include <stan/math/prim/err/make_iter_name.hpp>
+#include <stan/math/prim/err/throw_domain_error.hpp>
+#include <stan/math/prim/fun/to_ref.hpp>
+#include <stan/math/prim/fun/value_of_rec.hpp>
+#include <sstream>
+#include <string>
+
+namespace stan {
+namespace math {
+
+/**
+ * Throw an exception if the specified vector does not sum to 0.
+ * This function tests that the sum is within the tolerance specified by
+ * `CONSTRAINT_TOLERANCE`.
+ * This function only accepts Eigen vectors, statically
+ * typed vectors, not general matrices with 1 column.
+ * @tparam T A type inheriting from `Eigen::EigenBase`
+ * @param function Function name (for error messages)
+ * @param name Variable name (for error messages)
+ * @param theta Vector to test
+ * @throw `std::invalid_argument` if `theta` is a 0-vector
+ * @throw `std::domain_error` if the vector does not sum to zero
+ */
+template <typename T, require_matrix_t<T>* = nullptr>
+void check_sum_to_zero(const char* function, const char* name, const T& theta) {
+  using std::fabs;
+  // the size-zero case is technically a valid sum-to-zero vector,
+  // but it cannot be unconstrained to anything
+  check_nonzero_size(function, name, theta);
+  auto&& theta_ref = to_ref(value_of_rec(theta));
+  if (unlikely(!(fabs(theta_ref.sum()) <= CONSTRAINT_TOLERANCE))) {
+    [&]() STAN_COLD_PATH {
+      std::stringstream msg;
+      scalar_type_t<T> sum = theta_ref.sum();
+      msg << "does not sum to zero.";
+      msg.precision(10);
+      msg << " sum(" << name << ") = " << sum << ", but should be ";
+      std::string msg_str(msg.str());
+      throw_domain_error(function, name, 0.0, msg_str.c_str());
+    }();
+  }
+}
+
+/**
+ * Throw an exception if any vector in a standard vector does not sum to 0.
+ * This function tests that the sum is within the tolerance specified by
+ * `CONSTRAINT_TOLERANCE`.
+ * @tparam T A standard vector with inner type inheriting from
+ * `Eigen::EigenBase`
+ * @param function Function name (for error messages)
+ * @param name Variable name (for error messages)
+ * @param theta Vector to test.
+ * @throw `std::invalid_argument` if `theta` is a 0-vector
+ * @throw `std::domain_error` if the vector does not sum to zero
+ */
+template <typename T, require_std_vector_t<T>* = nullptr>
+void check_sum_to_zero(const char* function, const char* name, const T& theta) {
+  for (size_t i = 0; i < theta.size(); ++i) {
+    check_sum_to_zero(function, internal::make_iter_name(name, i).c_str(),
+                      theta[i]);
+  }
+}
+
+}  // namespace math
+}  // namespace stan
+#endif
diff --git a/stan/math/rev/constraint.hpp b/stan/math/rev/constraint.hpp
index 3ab4a6f72c7..ac036a31f52 100644
--- a/stan/math/rev/constraint.hpp
+++ b/stan/math/rev/constraint.hpp
@@ -17,6 +17,7 @@
 #include <stan/math/rev/constraint/simplex_constrain.hpp>
 #include <stan/math/rev/constraint/stochastic_column_constrain.hpp>
 #include <stan/math/rev/constraint/stochastic_row_constrain.hpp>
+#include <stan/math/rev/constraint/sum_to_zero_constrain.hpp>
 #include <stan/math/rev/constraint/unit_vector_constrain.hpp>
 #include <stan/math/rev/constraint/ub_constrain.hpp>
 
diff --git a/stan/math/rev/constraint/sum_to_zero_constrain.hpp b/stan/math/rev/constraint/sum_to_zero_constrain.hpp
new file mode 100644
index 00000000000..38248fe1413
--- /dev/null
+++ b/stan/math/rev/constraint/sum_to_zero_constrain.hpp
@@ -0,0 +1,104 @@
+#ifndef STAN_MATH_REV_CONSTRAINT_SUM_TO_ZERO_CONSTRAIN_HPP
+#define STAN_MATH_REV_CONSTRAINT_SUM_TO_ZERO_CONSTRAIN_HPP
+
+#include <stan/math/rev/meta.hpp>
+#include <stan/math/rev/core/reverse_pass_callback.hpp>
+#include <stan/math/rev/core/arena_matrix.hpp>
+#include <stan/math/prim/fun/Eigen.hpp>
+#include <stan/math/prim/fun/sqrt.hpp>
+#include <stan/math/prim/constraint/sum_to_zero_constrain.hpp>
+#include <cmath>
+#include <tuple>
+#include <vector>
+
+namespace stan {
+namespace math {
+
+/**
+ * Return a vector with sum zero corresponding to the specified
+ * free vector.
+ *
+ * The sum-to-zero transform is defined using a modified version of
+ * the inverse of the isometric log ratio transform (ILR).
+ * See:
+ * Egozcue, Juan Jose; Pawlowsky-Glahn, Vera; Mateu-Figueras, Gloria;
+ * Barcelo-Vidal, Carles (2003), "Isometric logratio transformations for
+ * compositional data analysis", Mathematical Geology, 35 (3): 279–300,
+ * doi:10.1023/A:1023818214614, S2CID 122844634
+ *
+ * This implementation is closer to the description of the same using "pivot
+ * coordinates" in
+ * Filzmoser, P., Hron, K., Templ, M. (2018). Geometrical Properties of
+ * Compositional Data. In: Applied Compositional Data Analysis. Springer Series
+ * in Statistics. Springer, Cham. https://doi.org/10.1007/978-3-319-96422-5_3
+ *
+ * This is a linear transform, with no Jacobian.
+ *
+ * @tparam T type of the vector
+ * @param y Free vector input of dimensionality K - 1.
+ * @return Zero-sum vector of dimensionality K.
+ */
+template <typename T, require_rev_col_vector_t<T>* = nullptr>
+inline auto sum_to_zero_constrain(T&& y) {
+  using ret_type = plain_type_t<T>;
+  if (unlikely(y.size() == 0)) {
+    return arena_t<ret_type>(Eigen::VectorXd{{0}});
+  }
+  auto arena_y = to_arena(std::forward<T>(y));
+  arena_t<ret_type> arena_z = sum_to_zero_constrain(arena_y.val());
+
+  reverse_pass_callback([arena_y, arena_z]() mutable {
+    const auto N = arena_y.size();
+
+    double sum_u_adj = 0;
+    for (int i = 0; i < N; ++i) {
+      double n = static_cast<double>(i + 1);
+
+      // adjoint of the reverse cumulative sum computed in the forward mode
+      sum_u_adj += arena_z.adj()(i);
+
+      // adjoint of the offset subtraction
+      double v_adj = -arena_z.adj()(i + 1) * n;
+
+      double w_adj = v_adj + sum_u_adj;
+
+      arena_y.adj()(i) += w_adj / sqrt(n * (n + 1));
+    }
+  });
+
+  return arena_z;
+}
+
+/**
+ * Return a vector with sum zero corresponding to the specified
+ * free vector.
+ *
+ * The sum-to-zero transform is defined using a modified version of
+ * the inverse of the isometric log ratio transform (ILR).
+ * See:
+ * Egozcue, Juan Jose; Pawlowsky-Glahn, Vera; Mateu-Figueras, Gloria;
+ * Barcelo-Vidal, Carles (2003), "Isometric logratio transformations for
+ * compositional data analysis", Mathematical Geology, 35 (3): 279–300,
+ * doi:10.1023/A:1023818214614, S2CID 122844634
+ *
+ * This implementation is closer to the description of the same using "pivot
+ * coordinates" in
+ * Filzmoser, P., Hron, K., Templ, M. (2018). Geometrical Properties of
+ * Compositional Data. In: Applied Compositional Data Analysis. Springer Series
+ * in Statistics. Springer, Cham. https://doi.org/10.1007/978-3-319-96422-5_3
+ *
+ * This is a linear transform, with no Jacobian.
+ *
+ * @tparam Vec type of the vector
+ * @param y Free vector input of dimensionality K - 1.
+ * @param lp unused
+ * @return Zero-sum vector of dimensionality K.
+ */
+template <typename T, require_rev_col_vector_t<T>* = nullptr>
+inline auto sum_to_zero_constrain(T&& y, scalar_type_t<T>& lp) {
+  return sum_to_zero_constrain(std::forward<T>(y));
+}
+
+}  // namespace math
+}  // namespace stan
+#endif
diff --git a/test/unit/math/mix/constraint/sum_to_zero_constrain_test.cpp b/test/unit/math/mix/constraint/sum_to_zero_constrain_test.cpp
new file mode 100644
index 00000000000..a42db039fcd
--- /dev/null
+++ b/test/unit/math/mix/constraint/sum_to_zero_constrain_test.cpp
@@ -0,0 +1,58 @@
+#include <test/unit/math/test_ad.hpp>
+
+namespace sum_to_zero_constrain_test {
+template <typename T>
+T g1(const T& x) {
+  stan::scalar_type_t<T> lp = 0;
+  return stan::math::sum_to_zero_constrain<false>(x, lp);
+}
+template <typename T>
+T g2(const T& x) {
+  stan::scalar_type_t<T> lp = 0;
+  return stan::math::sum_to_zero_constrain<true>(x, lp);
+}
+template <typename T>
+typename stan::scalar_type<T>::type g3(const T& x) {
+  stan::scalar_type_t<T> lp = 0;
+  stan::math::sum_to_zero_constrain<true>(x, lp);
+  return lp;
+}
+
+template <typename T>
+void expect_sum_to_zero_transform(const T& x) {
+  auto f1 = [](const auto& x) { return g1(x); };
+  auto f2 = [](const auto& x) { return g2(x); };
+  auto f3 = [](const auto& x) { return g3(x); };
+  stan::test::expect_ad(f1, x);
+  stan::test::expect_ad_matvar(f1, x);
+  stan::test::expect_ad(f2, x);
+  stan::test::expect_ad_matvar(f2, x);
+  stan::test::expect_ad(f3, x);
+  stan::test::expect_ad_matvar(f3, x);
+}
+}  // namespace sum_to_zero_constrain_test
+
+TEST(MathMixMatFun, sum_to_zeroTransform) {
+  Eigen::VectorXd v0(0);
+  sum_to_zero_constrain_test::expect_sum_to_zero_transform(v0);
+
+  Eigen::VectorXd v1(1);
+  v1 << 1;
+  sum_to_zero_constrain_test::expect_sum_to_zero_transform(v1);
+
+  Eigen::VectorXd v2(2);
+  v2 << 3, -1;
+  sum_to_zero_constrain_test::expect_sum_to_zero_transform(v2);
+
+  Eigen::VectorXd v3(3);
+  v3 << 2, 3, -1;
+  sum_to_zero_constrain_test::expect_sum_to_zero_transform(v3);
+
+  Eigen::VectorXd v4(4);
+  v4 << 2, -1, 0, -1.1;
+  sum_to_zero_constrain_test::expect_sum_to_zero_transform(v4);
+
+  Eigen::VectorXd v5(5);
+  v5 << 1, -3, 2, 0, -1;
+  sum_to_zero_constrain_test::expect_sum_to_zero_transform(v5);
+}
diff --git a/test/unit/math/prim/constraint/sum_to_zero_transform_test.cpp b/test/unit/math/prim/constraint/sum_to_zero_transform_test.cpp
new file mode 100644
index 00000000000..3c09aa1b768
--- /dev/null
+++ b/test/unit/math/prim/constraint/sum_to_zero_transform_test.cpp
@@ -0,0 +1,64 @@
+#include <stan/math/prim.hpp>
+#include <test/unit/util.hpp>
+#include <gtest/gtest.h>
+
+TEST(prob_transform, sum_to_zero_rt0) {
+  using Eigen::Dynamic;
+  using Eigen::Matrix;
+  double lp = 0;
+  Matrix<double, Dynamic, 1> x(4);
+  x << 0.0, 0.0, 0.0, 0.0;
+  std::vector<Matrix<double, Dynamic, 1>> x_vec{x, x, x};
+  std::vector<Matrix<double, Dynamic, 1>> y_vec
+      = stan::math::sum_to_zero_constrain<false>(x_vec, lp);
+  EXPECT_NO_THROW(stan::math::check_sum_to_zero("checkSumToZero", "y", y_vec));
+
+  for (auto&& y_i : y_vec) {
+    EXPECT_MATRIX_FLOAT_EQ(Eigen::VectorXd::Zero(5), y_i);
+  }
+  std::vector<Matrix<double, Dynamic, 1>> xrt
+      = stan::math::sum_to_zero_free(y_vec);
+  EXPECT_EQ(x.size() + 1, y_vec[2].size());
+  for (auto&& x_i : xrt) {
+    EXPECT_EQ(x.size(), x_i.size());
+    for (int i = 0; i < x.size(); ++i) {
+      EXPECT_NEAR(x[i], x_i[i], 1E-10);
+    }
+  }
+}
+TEST(prob_transform, sum_to_zero_rt) {
+  using Eigen::Dynamic;
+  using Eigen::Matrix;
+  Matrix<double, Dynamic, 1> x(3);
+  x << 1.0, -1.0, 2.0;
+  Matrix<double, Dynamic, 1> y = stan::math::sum_to_zero_constrain(x);
+  EXPECT_NO_THROW(stan::math::check_sum_to_zero("checkSumToZero", "y", y));
+  Matrix<double, Dynamic, 1> xrt = stan::math::sum_to_zero_free(y);
+  EXPECT_EQ(x.size() + 1, y.size());
+  EXPECT_EQ(x.size(), xrt.size());
+  for (int i = 0; i < x.size(); ++i) {
+    EXPECT_FLOAT_EQ(x[i], xrt[i]);
+  }
+}
+TEST(prob_transform, sum_to_zero_match) {
+  using Eigen::Dynamic;
+  using Eigen::Matrix;
+  Matrix<double, Dynamic, 1> x(3);
+  x << 1.0, -1.0, 2.0;
+  double lp = 0;
+  Matrix<double, Dynamic, 1> y = stan::math::sum_to_zero_constrain(x);
+  Matrix<double, Dynamic, 1> y2 = stan::math::sum_to_zero_constrain(x, lp);
+
+  EXPECT_EQ(4, y.size());
+  EXPECT_EQ(4, y2.size());
+  for (int i = 0; i < x.size(); ++i)
+    EXPECT_FLOAT_EQ(y[i], y2[i]);
+}
+
+TEST(prob_transform, sum_to_zero_f_exception) {
+  using Eigen::Dynamic;
+  using Eigen::Matrix;
+  Matrix<double, Dynamic, 1> y(2);
+  y << 0.5, -0.55;
+  EXPECT_THROW(stan::math::sum_to_zero_free(y), std::domain_error);
+}
diff --git a/test/unit/math/prim/err/check_sum_to_zero_test.cpp b/test/unit/math/prim/err/check_sum_to_zero_test.cpp
new file mode 100644
index 00000000000..1e6ad716a81
--- /dev/null
+++ b/test/unit/math/prim/err/check_sum_to_zero_test.cpp
@@ -0,0 +1,92 @@
+#include <stan/math/prim.hpp>
+#include <gtest/gtest.h>
+#include <test/unit/util.hpp>
+#include <limits>
+#include <string>
+
+TEST(ErrorHandlingMatrix, checkSumToZero_edges) {
+  Eigen::Matrix<double, Eigen::Dynamic, 1> zero(0);
+
+  EXPECT_THROW(stan::math::check_sum_to_zero("checkSumToZero", "zero", zero),
+               std::invalid_argument);
+
+  Eigen::Matrix<double, Eigen::Dynamic, 1> y_vec(1);
+  y_vec << 0.0;
+
+  EXPECT_NO_THROW(stan::math::check_sum_to_zero("checkSumToZero", "y", y_vec));
+
+  y_vec[0] = 0.1;
+  EXPECT_THROW(stan::math::check_sum_to_zero("checkSumToZero", "y", y_vec),
+               std::domain_error);
+}
+
+TEST(ErrorHandlingMatrix, checkSumToZero) {
+  Eigen::Matrix<double, Eigen::Dynamic, 1> y_vec(2);
+  std::vector<Eigen::Matrix<double, Eigen::Dynamic, 1>> y{y_vec, y_vec, y_vec};
+  for (auto& y_i : y) {
+    y_i << 0.5, -0.5;
+  }
+
+  EXPECT_NO_THROW(stan::math::check_sum_to_zero("checkSumToZero", "y", y));
+
+  for (auto& y_i : y) {
+    y_i[0] = 0.55;
+  }
+  EXPECT_THROW(stan::math::check_sum_to_zero("checkSumToZero", "y", y),
+               std::domain_error);
+}
+
+TEST(ErrorHandlingMatrix, checkSumToZero_message_sum) {
+  std::string message;
+  Eigen::Matrix<double, Eigen::Dynamic, 1> y_vec(100);
+  y_vec.setZero();
+  std::vector<Eigen::Matrix<double, Eigen::Dynamic, 1>> y{y_vec, y_vec, y_vec};
+  for (auto& y_i : y) {
+    y_i[13] = 0.1;
+  }
+
+  try {
+    stan::math::check_sum_to_zero("checkSumToZero", "y", y);
+    FAIL() << "should have thrown";
+  } catch (std::domain_error& e) {
+    message = e.what();
+  } catch (...) {
+    FAIL() << "threw the wrong error";
+  }
+
+  EXPECT_TRUE(std::string::npos != message.find(" y[1] does not sum to zero"))
+      << message;
+
+  EXPECT_TRUE(std::string::npos != message.find("sum(y[1]) = 0.1")) << message;
+}
+
+TEST(ErrorHandlingMatrix, checkSumToZero_nan) {
+  Eigen::Matrix<double, Eigen::Dynamic, 1> y_vec(2);
+  std::vector<Eigen::Matrix<double, Eigen::Dynamic, 1>> y{y_vec, y_vec, y_vec};
+  constexpr double nan = std::numeric_limits<double>::quiet_NaN();
+  for (auto& y_i : y) {
+    y_i << nan, 0.5;
+  }
+
+  EXPECT_THROW(stan::math::check_sum_to_zero("checkSumToZero", "y", y),
+               std::domain_error);
+
+  for (auto& y_i : y) {
+    y_i[1] = 0.55;
+  }
+  EXPECT_THROW(stan::math::check_sum_to_zero("checkSumToZero", "y", y),
+               std::domain_error);
+
+  for (auto& y_i : y) {
+    y_i[0] = 0.5;
+    y_i[1] = nan;
+  }
+  EXPECT_THROW(stan::math::check_sum_to_zero("checkSumToZero", "y", y),
+               std::domain_error);
+
+  for (auto& y_i : y) {
+    y_i[0] = nan;
+  }
+  EXPECT_THROW(stan::math::check_sum_to_zero("checkSumToZero", "y", y),
+               std::domain_error);
+}