diff --git a/stan/math/fwd.hpp b/stan/math/fwd.hpp
index e27b1b39915..36fa9e0e794 100644
--- a/stan/math/fwd.hpp
+++ b/stan/math/fwd.hpp
@@ -5,9 +5,9 @@
 
 #include <stan/math/fwd/core.hpp>
 #include <stan/math/fwd/meta.hpp>
-#include <stan/math/prim.hpp>
-
 #include <stan/math/fwd/fun.hpp>
 #include <stan/math/fwd/functor.hpp>
 
+#include <stan/math/prim.hpp>
+
 #endif
diff --git a/stan/math/fwd/fun/sum.hpp b/stan/math/fwd/fun/sum.hpp
index d2d4dca6222..2555b4f0122 100644
--- a/stan/math/fwd/fun/sum.hpp
+++ b/stan/math/fwd/fun/sum.hpp
@@ -1,9 +1,10 @@
 #ifndef STAN_MATH_FWD_FUN_SUM_HPP
 #define STAN_MATH_FWD_FUN_SUM_HPP
 
-#include <stan/math/fwd/core.hpp>
+#include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/fun/Eigen.hpp>
 #include <stan/math/prim/fun/sum.hpp>
+#include <stan/math/fwd/core.hpp>
 #include <vector>
 
 namespace stan {
diff --git a/stan/math/mix.hpp b/stan/math/mix.hpp
index b2814caa769..d60496114fa 100644
--- a/stan/math/mix.hpp
+++ b/stan/math/mix.hpp
@@ -2,11 +2,19 @@
 #define STAN_MATH_MIX_HPP
 
 #include <stan/math/mix/meta.hpp>
-
-#include <stan/math/prim.hpp>
-#include <stan/math/fwd.hpp>
-#include <stan/math/rev.hpp>
 #include <stan/math/mix/fun.hpp>
 #include <stan/math/mix/functor.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/fwd/core.hpp>
+#include <stan/math/fwd/meta.hpp>
+#include <stan/math/fwd/fun.hpp>
+#include <stan/math/fwd/functor.hpp>
+
+#include <stan/math/prim.hpp>
+
 #endif
diff --git a/stan/math/prim/err/elementwise_check.hpp b/stan/math/prim/err/elementwise_check.hpp
index effa475b83a..c6aaa127f93 100644
--- a/stan/math/prim/err/elementwise_check.hpp
+++ b/stan/math/prim/err/elementwise_check.hpp
@@ -142,8 +142,9 @@ inline void elementwise_check(const F& is_good, const char* function,
  */
 template <typename F, typename T, typename... Indexings,
           require_eigen_t<T>* = nullptr,
-          std::enable_if_t<static_cast<bool>(Eigen::internal::traits<T>::Flags&(
-              Eigen::LinearAccessBit | Eigen::DirectAccessBit))>* = nullptr>
+          std::enable_if_t<(Eigen::internal::traits<T>::Flags
+                            & Eigen::LinearAccessBit)
+                           || T::IsVectorAtCompileTime>* = nullptr>
 inline void elementwise_check(const F& is_good, const char* function,
                               const char* name, const T& x, const char* must_be,
                               const Indexings&... indexings) {
@@ -189,13 +190,13 @@ inline void elementwise_check(const F& is_good, const char* function,
  * @throws `std::domain_error` if `is_good` returns `false` for the value
  * of any element in `x`
  */
-template <
-    typename F, typename T, typename... Indexings,
-    require_eigen_t<T>* = nullptr,
-    std::enable_if_t<!(Eigen::internal::traits<T>::Flags
-                       & (Eigen::LinearAccessBit | Eigen::DirectAccessBit))
-                     && !(Eigen::internal::traits<T>::Flags
-                          & Eigen::RowMajorBit)>* = nullptr>
+template <typename F, typename T, typename... Indexings,
+          require_eigen_t<T>* = nullptr,
+          std::enable_if_t<!(Eigen::internal::traits<T>::Flags
+                             & Eigen::LinearAccessBit)
+                           && !T::IsVectorAtCompileTime
+                           && !(Eigen::internal::traits<T>::Flags
+                                & Eigen::RowMajorBit)>* = nullptr>
 inline void elementwise_check(const F& is_good, const char* function,
                               const char* name, const T& x, const char* must_be,
                               const Indexings&... indexings) {
@@ -230,13 +231,13 @@ inline void elementwise_check(const F& is_good, const char* function,
  * @throws `std::domain_error` if `is_good` returns `false` for the value
  * of any element in `x`
  */
-template <
-    typename F, typename T, typename... Indexings,
-    require_eigen_t<T>* = nullptr,
-    std::enable_if_t<!(Eigen::internal::traits<T>::Flags
-                       & (Eigen::LinearAccessBit | Eigen::DirectAccessBit))
-                     && static_cast<bool>(Eigen::internal::traits<T>::Flags
-                                          & Eigen::RowMajorBit)>* = nullptr>
+template <typename F, typename T, typename... Indexings,
+          require_eigen_t<T>* = nullptr,
+          std::enable_if_t<
+              !(Eigen::internal::traits<T>::Flags & Eigen::LinearAccessBit)
+              && !T::IsVectorAtCompileTime
+              && static_cast<bool>(Eigen::internal::traits<T>::Flags
+                                   & Eigen::RowMajorBit)>* = nullptr>
 inline void elementwise_check(const F& is_good, const char* function,
                               const char* name, const T& x, const char* must_be,
                               const Indexings&... indexings) {
diff --git a/stan/math/prim/fun/add.hpp b/stan/math/prim/fun/add.hpp
index 515e4d9ba90..763c51da1e9 100644
--- a/stan/math/prim/fun/add.hpp
+++ b/stan/math/prim/fun/add.hpp
@@ -42,7 +42,7 @@ template <typename Mat1, typename Mat2,
           require_all_not_st_var<Mat1, Mat2>* = nullptr>
 inline auto add(const Mat1& m1, const Mat2& m2) {
   check_matching_dims("add", "m1", m1, "m2", m2);
-  return (m1 + m2).eval();
+  return m1 + m2;
 }
 
 /**
@@ -58,7 +58,7 @@ template <typename Mat, typename Scal, require_eigen_t<Mat>* = nullptr,
           require_stan_scalar_t<Scal>* = nullptr,
           require_all_not_st_var<Mat, Scal>* = nullptr>
 inline auto add(const Mat& m, const Scal c) {
-  return (m.array() + c).matrix().eval();
+  return (m.array() + c).matrix();
 }
 
 /**
@@ -74,7 +74,7 @@ template <typename Scal, typename Mat, require_stan_scalar_t<Scal>* = nullptr,
           require_eigen_t<Mat>* = nullptr,
           require_all_not_st_var<Scal, Mat>* = nullptr>
 inline auto add(const Scal c, const Mat& m) {
-  return (c + m.array()).matrix().eval();
+  return (c + m.array()).matrix();
 }
 
 }  // namespace math
diff --git a/stan/math/prim/fun/block.hpp b/stan/math/prim/fun/block.hpp
index 9993412ace6..675a474fdc3 100644
--- a/stan/math/prim/fun/block.hpp
+++ b/stan/math/prim/fun/block.hpp
@@ -19,8 +19,7 @@ namespace math {
  * @throw std::out_of_range if either index is out of range.
  */
 template <typename T, require_eigen_t<T>* = nullptr>
-inline Eigen::Matrix<value_type_t<T>, Eigen::Dynamic, Eigen::Dynamic> block(
-    const T& m, size_t i, size_t j, size_t nrows, size_t ncols) {
+inline auto block(const T& m, size_t i, size_t j, size_t nrows, size_t ncols) {
   check_row_index("block", "i", m, i);
   check_row_index("block", "i+nrows-1", m, i + nrows - 1);
   check_column_index("block", "j", m, j);
diff --git a/stan/math/prim/fun/col.hpp b/stan/math/prim/fun/col.hpp
index d4f0298eccf..ee83247e1f2 100644
--- a/stan/math/prim/fun/col.hpp
+++ b/stan/math/prim/fun/col.hpp
@@ -23,7 +23,7 @@ namespace math {
 template <typename T, typename = require_eigen_t<T>>
 inline auto col(const T& m, size_t j) {
   check_column_index("col", "j", m, j);
-  return m.col(j - 1).eval();
+  return m.col(j - 1);
 }
 
 }  // namespace math
diff --git a/stan/math/prim/fun/corr_matrix_free.hpp b/stan/math/prim/fun/corr_matrix_free.hpp
index b17e4cff800..3d001173634 100644
--- a/stan/math/prim/fun/corr_matrix_free.hpp
+++ b/stan/math/prim/fun/corr_matrix_free.hpp
@@ -40,16 +40,16 @@ Eigen::Matrix<value_type_t<T>, Eigen::Dynamic, 1> corr_matrix_free(const T& y) {
 
   Eigen::Index k = y.rows();
   Eigen::Index k_choose_2 = (k * (k - 1)) / 2;
-  Array<value_type_t<T>, Dynamic, 1> x(k_choose_2);
+  Eigen::Matrix<value_type_t<T>, Dynamic, 1> x(k_choose_2);
   Array<value_type_t<T>, Dynamic, 1> sds(k);
-  bool successful = factor_cov_matrix(y, x, sds);
+  bool successful = factor_cov_matrix(y, x.array(), sds);
   if (!successful) {
     throw_domain_error("corr_matrix_free", "factor_cov_matrix failed on y", y,
                        "");
   }
   check_bounded("corr_matrix_free", "log(sd)", sds, -CONSTRAINT_TOLERANCE,
                 CONSTRAINT_TOLERANCE);
-  return x.matrix();
+  return x;
 }
 
 }  // namespace math
diff --git a/stan/math/prim/fun/diag_post_multiply.hpp b/stan/math/prim/fun/diag_post_multiply.hpp
index bfdae96b278..84ac4411cec 100644
--- a/stan/math/prim/fun/diag_post_multiply.hpp
+++ b/stan/math/prim/fun/diag_post_multiply.hpp
@@ -12,7 +12,7 @@ template <typename T1, typename T2, require_eigen_t<T1>* = nullptr,
 auto diag_post_multiply(const T1& m1, const T2& m2) {
   check_size_match("diag_post_multiply", "m2.size()", m2.size(), "m1.cols()",
                    m1.cols());
-  return (m1 * m2.asDiagonal()).eval();
+  return m1 * m2.asDiagonal();
 }
 
 }  // namespace math
diff --git a/stan/math/prim/fun/diag_pre_multiply.hpp b/stan/math/prim/fun/diag_pre_multiply.hpp
index 814378339c2..3302b75c086 100644
--- a/stan/math/prim/fun/diag_pre_multiply.hpp
+++ b/stan/math/prim/fun/diag_pre_multiply.hpp
@@ -12,7 +12,7 @@ template <typename T1, typename T2, require_eigen_vector_t<T1>* = nullptr,
 auto diag_pre_multiply(const T1& m1, const T2& m2) {
   check_size_match("diag_pre_multiply", "m1.size()", m1.size(), "m2.rows()",
                    m2.rows());
-  return (m1.asDiagonal() * m2).eval();
+  return m1.asDiagonal() * m2;
 }
 
 }  // namespace math
diff --git a/stan/math/prim/fun/diagonal.hpp b/stan/math/prim/fun/diagonal.hpp
index 560f680bd08..4b54385ca68 100644
--- a/stan/math/prim/fun/diagonal.hpp
+++ b/stan/math/prim/fun/diagonal.hpp
@@ -17,7 +17,7 @@ namespace math {
  */
 template <typename T, typename = require_eigen_t<T>>
 inline auto diagonal(const T& m) {
-  return m.diagonal().eval();
+  return m.diagonal();
 }
 
 }  // namespace math
diff --git a/stan/math/prim/fun/divide.hpp b/stan/math/prim/fun/divide.hpp
index fd7e48de087..c4d32fad6df 100644
--- a/stan/math/prim/fun/divide.hpp
+++ b/stan/math/prim/fun/divide.hpp
@@ -19,7 +19,7 @@ namespace math {
  * @return Scalar divided by the scalar.
  */
 template <typename Scal1, typename Scal2,
-          typename = require_all_stan_scalar_t<Scal1, Scal2>>
+          require_all_stan_scalar_t<Scal1, Scal2>* = nullptr>
 inline return_type_t<Scal1, Scal2> divide(const Scal1& x, const Scal2& y) {
   return x / y;
 }
@@ -41,10 +41,10 @@ inline int divide(int x, int y) {
  * @return matrix divided by the scalar
  */
 template <typename Mat, typename Scal, typename = require_eigen_t<Mat>,
-          typename = require_stan_scalar_t<Scal>,
-          typename = require_all_not_var_t<scalar_type_t<Mat>, Scal>>
+          require_stan_scalar_t<Scal>* = nullptr,
+          require_all_not_var_t<scalar_type_t<Mat>, Scal>* = nullptr>
 inline auto divide(const Mat& m, Scal c) {
-  return (m / c).eval();
+  return m / c;
 }
 
 }  // namespace math
diff --git a/stan/math/prim/fun/elt_divide.hpp b/stan/math/prim/fun/elt_divide.hpp
index a8793286d7a..d0fedb64bdc 100644
--- a/stan/math/prim/fun/elt_divide.hpp
+++ b/stan/math/prim/fun/elt_divide.hpp
@@ -24,7 +24,7 @@ template <typename Mat1, typename Mat2,
           require_all_not_st_var<Mat1, Mat2>* = nullptr>
 auto elt_divide(const Mat1& m1, const Mat2& m2) {
   check_matching_dims("elt_divide", "m1", m1, "m2", m2);
-  return (m1.array() / m2.array()).matrix().eval();
+  return (m1.array() / m2.array()).matrix();
 }
 
 /**
@@ -38,8 +38,8 @@ auto elt_divide(const Mat1& m1, const Mat2& m2) {
  * @param s scalar
  * @return Elementwise division of a scalar by matrix.
  */
-template <typename Mat, typename Scal, typename = require_matrix_t<Mat>,
-          typename = require_stan_scalar_t<Scal>>
+template <typename Mat, typename Scal, require_matrix_t<Mat>* = nullptr,
+          require_stan_scalar_t<Scal>* = nullptr>
 auto elt_divide(const Mat& m, Scal s) {
   return divide(m, s);
 }
@@ -55,10 +55,10 @@ auto elt_divide(const Mat& m, Scal s) {
  * @param m matrix or expression
  * @return Elementwise division of a scalar by matrix.
  */
-template <typename Scal, typename Mat, typename = require_stan_scalar_t<Scal>,
-          typename = require_eigen_t<Mat>>
+template <typename Scal, typename Mat, require_stan_scalar_t<Scal>* = nullptr,
+          require_eigen_t<Mat>* = nullptr>
 auto elt_divide(Scal s, const Mat& m) {
-  return (s / m.array()).matrix().eval();
+  return (s / m.array()).matrix();
 }
 
 template <typename Scal1, typename Scal2,
diff --git a/stan/math/prim/fun/elt_multiply.hpp b/stan/math/prim/fun/elt_multiply.hpp
index 379429adb1a..060dc73d11f 100644
--- a/stan/math/prim/fun/elt_multiply.hpp
+++ b/stan/math/prim/fun/elt_multiply.hpp
@@ -24,7 +24,7 @@ template <typename Mat1, typename Mat2,
           require_all_not_st_var<Mat1, Mat2>* = nullptr>
 auto elt_multiply(const Mat1& m1, const Mat2& m2) {
   check_matching_dims("elt_multiply", "m1", m1, "m2", m2);
-  return m1.cwiseProduct(m2).eval();
+  return m1.cwiseProduct(m2);
 }
 
 /**
diff --git a/stan/math/prim/fun/head.hpp b/stan/math/prim/fun/head.hpp
index afea8187a0d..ac3afdd3be7 100644
--- a/stan/math/prim/fun/head.hpp
+++ b/stan/math/prim/fun/head.hpp
@@ -23,7 +23,7 @@ inline auto head(const T& v, size_t n) {
   if (n != 0) {
     check_vector_index("head", "n", v, n);
   }
-  return v.head(n).eval();
+  return v.head(n);
 }
 
 /**
diff --git a/stan/math/prim/fun/identity_matrix.hpp b/stan/math/prim/fun/identity_matrix.hpp
index b6b0542122c..eeb4b815916 100644
--- a/stan/math/prim/fun/identity_matrix.hpp
+++ b/stan/math/prim/fun/identity_matrix.hpp
@@ -14,7 +14,7 @@ namespace math {
  * @return An identity matrix of size K.
  * @throw std::domain_error if K is negative.
  */
-inline Eigen::MatrixXd identity_matrix(int K) {
+inline auto identity_matrix(int K) {
   check_nonnegative("identity_matrix", "size", K);
   return Eigen::MatrixXd::Identity(K, K);
 }
diff --git a/stan/math/prim/fun/inv_square.hpp b/stan/math/prim/fun/inv_square.hpp
index 15e40cf4427..53911218721 100644
--- a/stan/math/prim/fun/inv_square.hpp
+++ b/stan/math/prim/fun/inv_square.hpp
@@ -4,8 +4,6 @@
 #include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/fun/inv.hpp>
 #include <stan/math/prim/fun/square.hpp>
-#include <stan/math/prim/fun/inv_square.hpp>
-#include <stan/math/prim/functor/apply_scalar_unary.hpp>
 #include <stan/math/prim/functor/apply_vector_unary.hpp>
 
 namespace stan {
diff --git a/stan/math/prim/fun/linspaced_row_vector.hpp b/stan/math/prim/fun/linspaced_row_vector.hpp
index f1faf6398b6..d215c53c3b4 100644
--- a/stan/math/prim/fun/linspaced_row_vector.hpp
+++ b/stan/math/prim/fun/linspaced_row_vector.hpp
@@ -22,7 +22,7 @@ namespace math {
  * @throw std::domain_error if K is negative, if low is nan or infinite,
  * if high is nan or infinite, or if high is less than low.
  */
-inline Eigen::RowVectorXd linspaced_row_vector(int K, double low, double high) {
+inline auto linspaced_row_vector(int K, double low, double high) {
   static const char* function = "linspaced_row_vector";
   check_nonnegative(function, "size", K);
   check_finite(function, "low", low);
diff --git a/stan/math/prim/fun/linspaced_vector.hpp b/stan/math/prim/fun/linspaced_vector.hpp
index 96f9b0c6226..e0ab8a42dcf 100644
--- a/stan/math/prim/fun/linspaced_vector.hpp
+++ b/stan/math/prim/fun/linspaced_vector.hpp
@@ -22,7 +22,7 @@ namespace math {
  * @throw std::domain_error if K is negative, if low is nan or infinite,
  * if high is nan or infinite, or if high is less than low.
  */
-inline Eigen::VectorXd linspaced_vector(int K, double low, double high) {
+inline auto linspaced_vector(int K, double low, double high) {
   static const char* function = "linspaced_vector";
   check_nonnegative(function, "size", K);
   check_finite(function, "low", low);
diff --git a/stan/math/prim/fun/log1p.hpp b/stan/math/prim/fun/log1p.hpp
index 47b315a2bd7..94557fdd5a5 100644
--- a/stan/math/prim/fun/log1p.hpp
+++ b/stan/math/prim/fun/log1p.hpp
@@ -5,7 +5,6 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/Eigen.hpp>
 #include <stan/math/prim/fun/is_nan.hpp>
-#include <stan/math/prim/fun/log1p.hpp>
 #include <stan/math/prim/functor/apply_scalar_unary.hpp>
 #include <cmath>
 
diff --git a/stan/math/prim/fun/log_softmax.hpp b/stan/math/prim/fun/log_softmax.hpp
index f90ccee11b9..c863f23e612 100644
--- a/stan/math/prim/fun/log_softmax.hpp
+++ b/stan/math/prim/fun/log_softmax.hpp
@@ -5,6 +5,7 @@
 #include <stan/math/prim/err.hpp>
 #include <stan/math/prim/fun/Eigen.hpp>
 #include <stan/math/prim/fun/log_sum_exp.hpp>
+#include <stan/math/prim/fun/to_ref.hpp>
 #include <stan/math/prim/functor/apply_vector_unary.hpp>
 
 namespace stan {
@@ -37,18 +38,16 @@ namespace math {
  * @tparam Container type of input vector to transform
  * @param[in] x vector to transform
  * @return log unit simplex result of the softmax transform of the vector.
- *
- * Note: The return must be evaluated otherwise the Ref object falls out
- * of scope
  */
-template <typename Container,
-          require_arithmetic_t<scalar_type_t<Container>>* = nullptr>
+template <typename Container, require_st_arithmetic<Container>* = nullptr>
 inline auto log_softmax(const Container& x) {
-  return apply_vector_unary<Container>::apply(x, [](const auto& v) {
-    const Eigen::Ref<const plain_type_t<decltype(v)>>& v_ref = v;
-    check_nonzero_size("log_softmax", "v", v_ref);
-    return (v_ref.array() - log_sum_exp(v_ref)).eval();
-  });
+  check_nonzero_size("log_softmax", "v", x);
+  return make_holder(
+      [](const auto& a) {
+        return apply_vector_unary<ref_type_t<Container>>::apply(
+            a, [](const auto& v) { return v.array() - log_sum_exp(v); });
+      },
+      to_ref(x));
 }
 
 }  // namespace math
diff --git a/stan/math/prim/fun/logit.hpp b/stan/math/prim/fun/logit.hpp
index 2dc14698e80..8f71079217a 100644
--- a/stan/math/prim/fun/logit.hpp
+++ b/stan/math/prim/fun/logit.hpp
@@ -3,6 +3,7 @@
 
 #include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/fun/log.hpp>
+#include <stan/math/prim/fun/to_ref.hpp>
 #include <stan/math/prim/functor/apply_scalar_unary.hpp>
 #include <stan/math/prim/functor/apply_vector_unary.hpp>
 #include <cmath>
@@ -102,10 +103,13 @@ inline auto logit(const Container& x) {
 template <typename Container,
           require_container_st<std::is_arithmetic, Container>* = nullptr>
 inline auto logit(const Container& x) {
-  return apply_vector_unary<Container>::apply(x, [](const auto& v) {
-    const Eigen::Ref<const plain_type_t<decltype(v)>>& v_ref = v;
-    return (v_ref.array() / (1 - v_ref.array())).log().eval();
-  });
+  return make_holder(
+      [](const auto& v_ref) {
+        return apply_vector_unary<ref_type_t<Container>>::apply(
+            v_ref,
+            [](const auto& v) { return (v.array() / (1 - v.array())).log(); });
+      },
+      to_ref(x));
 }
 
 }  // namespace math
diff --git a/stan/math/prim/fun/minus.hpp b/stan/math/prim/fun/minus.hpp
index 13fe7d090ed..a7a8c2bf1b5 100644
--- a/stan/math/prim/fun/minus.hpp
+++ b/stan/math/prim/fun/minus.hpp
@@ -14,7 +14,7 @@ namespace math {
  * @return Negation of subtrahend.
  */
 template <typename T>
-inline plain_type_t<T> minus(const T& x) {
+inline auto minus(const T& x) {
   return -x;
 }
 
diff --git a/stan/math/prim/fun/multiply.hpp b/stan/math/prim/fun/multiply.hpp
index 59d2b37d8d5..d7f54687dfd 100644
--- a/stan/math/prim/fun/multiply.hpp
+++ b/stan/math/prim/fun/multiply.hpp
@@ -27,7 +27,7 @@ template <typename Mat, typename Scal, require_stan_scalar_t<Scal>* = nullptr,
           require_eigen_t<Mat>* = nullptr,
           require_all_not_st_var<Scal, Mat>* = nullptr>
 inline auto multiply(const Mat& m, Scal c) {
-  return (c * m).eval();
+  return c * m;
 }
 
 /**
@@ -44,7 +44,7 @@ template <typename Scal, typename Mat, require_stan_scalar_t<Scal>* = nullptr,
           require_eigen_t<Mat>* = nullptr,
           require_all_not_st_var<Scal, Mat>* = nullptr>
 inline auto multiply(Scal c, const Mat& m) {
-  return (c * m).eval();
+  return c * m;
 }
 
 /**
@@ -69,7 +69,7 @@ template <typename Mat1, typename Mat2,
 inline auto multiply(const Mat1& m1, const Mat2& m2) {
   check_size_match("multiply", "Columns of m1", m1.cols(), "Rows of m2",
                    m2.rows());
-  return (m1 * m2).eval();
+  return m1 * m2;
 }
 
 /**
diff --git a/stan/math/prim/fun/ones_row_vector.hpp b/stan/math/prim/fun/ones_row_vector.hpp
index 4bafd822f9e..d85b4364a38 100644
--- a/stan/math/prim/fun/ones_row_vector.hpp
+++ b/stan/math/prim/fun/ones_row_vector.hpp
@@ -14,7 +14,7 @@ namespace math {
  * @return A row vector of size K with all elements initialised to 1.
  * @throw std::domain_error if K is negative.
  */
-inline Eigen::RowVectorXd ones_row_vector(int K) {
+inline auto ones_row_vector(int K) {
   check_nonnegative("ones_row_vector", "size", K);
   return Eigen::RowVectorXd::Constant(K, 1);
 }
diff --git a/stan/math/prim/fun/ones_vector.hpp b/stan/math/prim/fun/ones_vector.hpp
index e6bf6b520ab..146c8e26f83 100644
--- a/stan/math/prim/fun/ones_vector.hpp
+++ b/stan/math/prim/fun/ones_vector.hpp
@@ -14,7 +14,7 @@ namespace math {
  * @return A vector of size K with all elements initialised to 1.
  * @throw std::domain_error if K is negative.
  */
-inline Eigen::VectorXd ones_vector(int K) {
+inline auto ones_vector(int K) {
   check_nonnegative("ones_vector", "size", K);
   return Eigen::VectorXd::Constant(K, 1);
 }
diff --git a/stan/math/prim/fun/quad_form.hpp b/stan/math/prim/fun/quad_form.hpp
index 1177fd66c92..f554d235eb2 100644
--- a/stan/math/prim/fun/quad_form.hpp
+++ b/stan/math/prim/fun/quad_form.hpp
@@ -28,13 +28,11 @@ template <typename EigMat1, typename EigMat2,
           require_not_eigen_col_vector_t<EigMat2>* = nullptr,
           require_vt_same<EigMat1, EigMat2>* = nullptr,
           require_all_vt_arithmetic<EigMat1, EigMat2>* = nullptr>
-inline Eigen::Matrix<value_type_t<EigMat2>, EigMat2::ColsAtCompileTime,
-                     EigMat2::ColsAtCompileTime>
-quad_form(const EigMat1& A, const EigMat2& B) {
+inline auto quad_form(const EigMat1& A, const EigMat2& B) {
   check_square("quad_form", "A", A);
   check_multiplicable("quad_form", "A", A, "B", B);
-  const auto& B_ref = to_ref(B);
-  return B_ref.transpose() * A * B_ref;
+  return make_holder([&A](const auto& b) { return b.transpose() * A * b; },
+                     to_ref(B));
 }
 
 /**
diff --git a/stan/math/prim/fun/quad_form_diag.hpp b/stan/math/prim/fun/quad_form_diag.hpp
index 7b243bb4e9b..093bef2275c 100644
--- a/stan/math/prim/fun/quad_form_diag.hpp
+++ b/stan/math/prim/fun/quad_form_diag.hpp
@@ -10,14 +10,13 @@ namespace math {
 
 template <typename EigMat, typename EigVec, require_eigen_t<EigMat>* = nullptr,
           require_eigen_vector_t<EigVec>* = nullptr>
-inline Eigen::Matrix<return_type_t<EigMat, EigVec>, Eigen::Dynamic,
-                     Eigen::Dynamic>
-quad_form_diag(const EigMat& mat, const EigVec& vec) {
+inline auto quad_form_diag(const EigMat& mat, const EigVec& vec) {
   check_square("quad_form_diag", "mat", mat);
   check_size_match("quad_form_diag", "rows of mat", mat.rows(), "size of vec",
                    vec.size());
-  const auto& vec_ref = to_ref(vec);
-  return vec_ref.asDiagonal() * mat * vec_ref.asDiagonal();
+  return make_holder(
+      [&mat](const auto& v) { return v.asDiagonal() * mat * v.asDiagonal(); },
+      to_ref(vec));
 }
 
 }  // namespace math
diff --git a/stan/math/prim/fun/quad_form_sym.hpp b/stan/math/prim/fun/quad_form_sym.hpp
index 7c367d5dda6..7b9359d82ae 100644
--- a/stan/math/prim/fun/quad_form_sym.hpp
+++ b/stan/math/prim/fun/quad_form_sym.hpp
@@ -32,8 +32,9 @@ inline plain_type_t<EigMat2> quad_form_sym(const EigMat1& A, const EigMat2& B) {
   const auto& A_ref = to_ref(A);
   const auto& B_ref = to_ref(B);
   check_symmetric("quad_form_sym", "A", A_ref);
-  plain_type_t<EigMat2> ret(B_ref.transpose() * A_ref * B_ref);
-  return value_type_t<EigMat2>(0.5) * (ret + ret.transpose());
+  return make_holder(
+      [](const auto& ret) { return 0.5 * (ret + ret.transpose()); },
+      (B_ref.transpose() * A_ref * B_ref).eval());
 }
 
 /**
diff --git a/stan/math/prim/fun/read_corr_L.hpp b/stan/math/prim/fun/read_corr_L.hpp
index 4d823de1d35..659bc872b0a 100644
--- a/stan/math/prim/fun/read_corr_L.hpp
+++ b/stan/math/prim/fun/read_corr_L.hpp
@@ -51,7 +51,7 @@ Eigen::Matrix<value_type_t<T>, Eigen::Dynamic, Eigen::Dynamic> read_corr_L(
   Eigen::Array<T_scalar, Eigen::Dynamic, 1> acc(K - 1);
   acc.setOnes();
   // Cholesky factor of correlation matrix
-  Eigen::Array<T_scalar, Eigen::Dynamic, Eigen::Dynamic> L(K, K);
+  Eigen::Matrix<T_scalar, Eigen::Dynamic, Eigen::Dynamic> L(K, K);
   L.setZero();
 
   size_t position = 0;
@@ -69,7 +69,7 @@ Eigen::Matrix<value_type_t<T>, Eigen::Dynamic, Eigen::Dynamic> read_corr_L(
     acc.tail(pull) *= T_scalar(1.0) - temp.square();
   }
   L(K - 1, K - 1) = sqrt(acc(K - 2));
-  return L.matrix();
+  return L;
 }
 
 /**
diff --git a/stan/math/prim/fun/read_cov_L.hpp b/stan/math/prim/fun/read_cov_L.hpp
index 17a6e19818d..6dce25df5d4 100644
--- a/stan/math/prim/fun/read_cov_L.hpp
+++ b/stan/math/prim/fun/read_cov_L.hpp
@@ -32,7 +32,9 @@ Eigen::Matrix<value_type<T_CPCs>, Eigen::Dynamic, Eigen::Dynamic> read_cov_L(
   size_t K = sds.rows();
   // adjust due to transformation from correlations to covariances
   log_prob += (sum(log(sds)) + LOG_TWO) * K;
-  return sds.matrix().asDiagonal() * read_corr_L(CPCs, K, log_prob);
+  return make_holder(
+      [&sds](const auto& b) { return sds.matrix().asDiagonal() * b; },
+      read_corr_L(CPCs, K, log_prob));
 }
 
 }  // namespace math
diff --git a/stan/math/prim/fun/rep_matrix.hpp b/stan/math/prim/fun/rep_matrix.hpp
index 4e3f20a9aa7..a3a3118fb50 100644
--- a/stan/math/prim/fun/rep_matrix.hpp
+++ b/stan/math/prim/fun/rep_matrix.hpp
@@ -9,8 +9,7 @@ namespace stan {
 namespace math {
 
 template <typename T, require_stan_scalar_t<T>* = nullptr>
-inline Eigen::Matrix<return_type_t<T>, Eigen::Dynamic, Eigen::Dynamic>
-rep_matrix(const T& x, int m, int n) {
+inline auto rep_matrix(const T& x, int m, int n) {
   check_nonnegative("rep_matrix", "rows", m);
   check_nonnegative("rep_matrix", "cols", n);
   return Eigen::Matrix<return_type_t<T>, Eigen::Dynamic,
@@ -18,15 +17,13 @@ rep_matrix(const T& x, int m, int n) {
 }
 
 template <typename ColVec, require_eigen_col_vector_t<ColVec>* = nullptr>
-inline Eigen::Matrix<value_type_t<ColVec>, Eigen::Dynamic, Eigen::Dynamic>
-rep_matrix(const ColVec& v, int n) {
+inline auto rep_matrix(const ColVec& v, int n) {
   check_nonnegative("rep_matrix", "rows", n);
   return v.replicate(1, n);
 }
 
 template <typename RowVec, require_eigen_row_vector_t<RowVec>* = nullptr>
-inline Eigen::Matrix<value_type_t<RowVec>, Eigen::Dynamic, Eigen::Dynamic>
-rep_matrix(const RowVec& rv, int m) {
+inline auto rep_matrix(const RowVec& rv, int m) {
   check_nonnegative("rep_matrix", "cols", m);
   return rv.replicate(m, 1);
 }
diff --git a/stan/math/prim/fun/rep_row_vector.hpp b/stan/math/prim/fun/rep_row_vector.hpp
index 6b4a95887d1..b79efdce89c 100644
--- a/stan/math/prim/fun/rep_row_vector.hpp
+++ b/stan/math/prim/fun/rep_row_vector.hpp
@@ -8,8 +8,7 @@ namespace stan {
 namespace math {
 
 template <typename T>
-inline Eigen::Matrix<return_type_t<T>, 1, Eigen::Dynamic> rep_row_vector(
-    const T& x, int m) {
+inline auto rep_row_vector(const T& x, int m) {
   check_nonnegative("rep_row_vector", "m", m);
   return Eigen::Matrix<return_type_t<T>, 1, Eigen::Dynamic>::Constant(m, x);
 }
diff --git a/stan/math/prim/fun/rep_vector.hpp b/stan/math/prim/fun/rep_vector.hpp
index 9668228a808..d0e883b07af 100644
--- a/stan/math/prim/fun/rep_vector.hpp
+++ b/stan/math/prim/fun/rep_vector.hpp
@@ -8,8 +8,7 @@ namespace stan {
 namespace math {
 
 template <typename T>
-inline Eigen::Matrix<return_type_t<T>, Eigen::Dynamic, 1> rep_vector(const T& x,
-                                                                     int n) {
+inline auto rep_vector(const T& x, int n) {
   check_nonnegative("rep_vector", "n", n);
   return Eigen::Matrix<return_type_t<T>, Eigen::Dynamic, 1>::Constant(n, x);
 }
diff --git a/stan/math/prim/fun/reverse.hpp b/stan/math/prim/fun/reverse.hpp
index 0fe1cbf5980..f8cf335eec7 100644
--- a/stan/math/prim/fun/reverse.hpp
+++ b/stan/math/prim/fun/reverse.hpp
@@ -33,7 +33,7 @@ inline std::vector<T> reverse(const std::vector<T>& x) {
  */
 template <typename T, typename = require_eigen_vector_t<T>>
 inline auto reverse(const T& x) {
-  return x.reverse().eval();
+  return x.reverse();
 }
 
 }  // namespace math
diff --git a/stan/math/prim/fun/row.hpp b/stan/math/prim/fun/row.hpp
index b9a142277b2..e51730dc90b 100644
--- a/stan/math/prim/fun/row.hpp
+++ b/stan/math/prim/fun/row.hpp
@@ -20,11 +20,11 @@ namespace math {
  * @return Specified row of the matrix.
  * @throw std::out_of_range if i is out of range.
  */
-template <typename T, typename = require_eigen_t<T>>
+template <typename T, require_eigen_t<T>* = nullptr>
 inline auto row(const T& m, size_t i) {
   check_row_index("row", "i", m, i);
 
-  return m.row(i - 1).eval();
+  return m.row(i - 1);
 }
 
 }  // namespace math
diff --git a/stan/math/prim/fun/segment.hpp b/stan/math/prim/fun/segment.hpp
index ec4b9aa9a70..109e37c9e8f 100644
--- a/stan/math/prim/fun/segment.hpp
+++ b/stan/math/prim/fun/segment.hpp
@@ -15,7 +15,7 @@ namespace math {
  * @tparam T type of the vector
  */
 template <typename EigVec, require_eigen_vector_t<EigVec>* = nullptr>
-inline plain_type_t<EigVec> segment(const EigVec& v, size_t i, size_t n) {
+inline auto segment(const EigVec& v, size_t i, size_t n) {
   check_greater("segment", "n", i, 0.0);
   check_less_or_equal("segment", "n", i, static_cast<size_t>(v.size()));
   if (n != 0) {
diff --git a/stan/math/prim/fun/softmax.hpp b/stan/math/prim/fun/softmax.hpp
index 88fbc3f2731..7a389559f99 100644
--- a/stan/math/prim/fun/softmax.hpp
+++ b/stan/math/prim/fun/softmax.hpp
@@ -50,8 +50,9 @@ inline plain_type_t<ColVec> softmax(const ColVec& v) {
     return v;
   }
   const auto& v_ref = to_ref(v);
-  plain_type_t<ColVec> theta = (v_ref.array() - v_ref.maxCoeff()).exp();
-  return theta.array() / theta.sum();
+  return make_holder(
+      [](const auto& theta) { return theta.array() / theta.sum(); },
+      (v_ref.array() - v_ref.maxCoeff()).exp().eval());
 }
 
 }  // namespace math
diff --git a/stan/math/prim/fun/sub_col.hpp b/stan/math/prim/fun/sub_col.hpp
index ab9c061ac15..a4fb70984a7 100644
--- a/stan/math/prim/fun/sub_col.hpp
+++ b/stan/math/prim/fun/sub_col.hpp
@@ -24,7 +24,7 @@ inline auto sub_col(const T& m, size_t i, size_t j, size_t nrows) {
     check_row_index("sub_col", "i+nrows-1", m, i + nrows - 1);
   }
   check_column_index("sub_col", "j", m, j);
-  return m.col(j - 1).segment(i - 1, nrows).eval();
+  return m.col(j - 1).segment(i - 1, nrows);
 }
 
 }  // namespace math
diff --git a/stan/math/prim/fun/sub_row.hpp b/stan/math/prim/fun/sub_row.hpp
index 15274656362..18b078cb015 100644
--- a/stan/math/prim/fun/sub_row.hpp
+++ b/stan/math/prim/fun/sub_row.hpp
@@ -24,7 +24,7 @@ inline auto sub_row(const T& m, size_t i, size_t j, size_t ncols) {
   if (ncols > 0) {
     check_column_index("sub_col", "j+ncols-1", m, j + ncols - 1);
   }
-  return m.row(i - 1).segment(j - 1, ncols).eval();
+  return m.row(i - 1).segment(j - 1, ncols);
 }
 
 }  // namespace math
diff --git a/stan/math/prim/fun/subtract.hpp b/stan/math/prim/fun/subtract.hpp
index 85ea3233158..33b369da108 100644
--- a/stan/math/prim/fun/subtract.hpp
+++ b/stan/math/prim/fun/subtract.hpp
@@ -42,7 +42,7 @@ template <typename Mat1, typename Mat2,
           require_all_not_st_var<Mat1, Mat2>* = nullptr>
 inline auto subtract(const Mat1& m1, const Mat2& m2) {
   check_matching_dims("subtract", "m1", m1, "m2", m2);
-  return (m1 - m2).eval();
+  return m1 - m2;
 }
 
 /**
@@ -59,7 +59,7 @@ template <typename Scal, typename Mat, typename = require_stan_scalar_t<Scal>,
           require_eigen_t<Mat>* = nullptr,
           require_all_not_st_var<Mat, Scal>* = nullptr>
 inline auto subtract(const Scal c, const Mat& m) {
-  return (c - m.array()).matrix().eval();
+  return (c - m.array()).matrix();
 }
 
 /**
@@ -72,11 +72,11 @@ inline auto subtract(const Scal c, const Mat& m) {
  * @param c Scalar.
  * @return The matrix minus the scalar.
  */
-template <typename Mat, typename Scal, typename = require_eigen_t<Mat>,
+template <typename Mat, typename Scal, require_eigen_t<Mat>* = nullptr,
           require_stan_scalar_t<Scal>* = nullptr,
           require_all_not_st_var<Scal, Mat>* = nullptr>
 inline auto subtract(const Mat& m, const Scal c) {
-  return (m.array() - c).matrix().eval();
+  return (m.array() - c).matrix();
 }
 
 }  // namespace math
diff --git a/stan/math/prim/fun/tail.hpp b/stan/math/prim/fun/tail.hpp
index bc7550083bc..75813acfddc 100644
--- a/stan/math/prim/fun/tail.hpp
+++ b/stan/math/prim/fun/tail.hpp
@@ -24,7 +24,7 @@ inline auto tail(const T& v, size_t n) {
   if (n != 0) {
     check_vector_index("tail", "n", v, n);
   }
-  return v.tail(n).eval();
+  return v.tail(n);
 }
 
 /**
diff --git a/stan/math/prim/fun/to_matrix.hpp b/stan/math/prim/fun/to_matrix.hpp
index a227c239d43..8f97574bc01 100644
--- a/stan/math/prim/fun/to_matrix.hpp
+++ b/stan/math/prim/fun/to_matrix.hpp
@@ -21,8 +21,7 @@ namespace math {
  * @return the matrix representation of the input
  */
 template <typename EigMat, require_eigen_t<EigMat>* = nullptr>
-inline Eigen::Matrix<value_type_t<EigMat>, Eigen::Dynamic, Eigen::Dynamic>
-to_matrix(EigMat&& x) {
+inline EigMat to_matrix(EigMat&& x) {
   return std::forward<EigMat>(x);
 }
 
@@ -114,8 +113,7 @@ to_matrix(EigMat&& x, int m, int n) {
  * if the sizes do not match
  */
 template <typename T>
-inline Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> to_matrix(
-    const std::vector<T>& x, int m, int n) {
+inline auto to_matrix(const std::vector<T>& x, int m, int n) {
   static const char* function = "to_matrix(array)";
   check_size_match(function, "rows * columns", m * n, "vector size", x.size());
   return Eigen::Map<const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>>(
diff --git a/stan/math/prim/fun/transpose.hpp b/stan/math/prim/fun/transpose.hpp
index 05a238d4884..70d860cb6cb 100644
--- a/stan/math/prim/fun/transpose.hpp
+++ b/stan/math/prim/fun/transpose.hpp
@@ -14,8 +14,7 @@ namespace math {
  * @return transposed matrix
  */
 template <typename T, typename = require_eigen_t<T>>
-Eigen::Matrix<value_type_t<T>, T::ColsAtCompileTime,
-              T::RowsAtCompileTime> inline transpose(const T& m) {
+auto inline transpose(const T& m) {
   return m.transpose();
 }
 
diff --git a/stan/math/prim/fun/uniform_simplex.hpp b/stan/math/prim/fun/uniform_simplex.hpp
index 7c473dbd90e..e2381160713 100644
--- a/stan/math/prim/fun/uniform_simplex.hpp
+++ b/stan/math/prim/fun/uniform_simplex.hpp
@@ -15,7 +15,7 @@ namespace math {
  * so that their sum is equal to 1.
  * @throw std::domain_error if K is not positive.
  */
-inline Eigen::VectorXd uniform_simplex(int K) {
+inline auto uniform_simplex(int K) {
   check_positive("uniform_simplex", "size", K);
   return Eigen::VectorXd::Constant(K, 1.0 / K);
 }
diff --git a/stan/math/prim/fun/unit_vector_constrain.hpp b/stan/math/prim/fun/unit_vector_constrain.hpp
index 718ddabb57c..d70ac629694 100644
--- a/stan/math/prim/fun/unit_vector_constrain.hpp
+++ b/stan/math/prim/fun/unit_vector_constrain.hpp
@@ -24,13 +24,16 @@ namespace math {
  */
 template <typename T, require_eigen_col_vector_t<T>* = nullptr,
           require_not_vt_autodiff<T>* = nullptr>
-inline plain_type_t<T> unit_vector_constrain(const T& y) {
+inline auto unit_vector_constrain(const T& y) {
   using std::sqrt;
   check_nonzero_size("unit_vector_constrain", "y", y);
-  const auto& y_ref = to_ref(y);
-  value_type_t<T> SN = dot_self(y_ref);
-  check_positive_finite("unit_vector_constrain", "norm", SN);
-  return (y_ref / sqrt(SN)).eval();
+  return make_holder(
+      [](const auto& y_ref) {
+        value_type_t<T> SN = dot_self(y_ref);
+        check_positive_finite("unit_vector_constrain", "norm", SN);
+        return y_ref / sqrt(SN);
+      },
+      to_ref(y));
 }
 
 /**
@@ -48,12 +51,15 @@ template <typename T1, typename T2, require_eigen_col_vector_t<T1>* = nullptr,
           require_all_not_vt_autodiff<T1, T2>* = nullptr>
 inline plain_type_t<T1> unit_vector_constrain(const T1& y, T2& lp) {
   using std::sqrt;
-  const auto& y_ref = to_ref(y);
-  check_nonzero_size("unit_vector_constrain", "y", y_ref);
-  value_type_t<T1> SN = dot_self(y_ref);
-  check_positive_finite("unit_vector_constrain", "norm", SN);
-  lp -= 0.5 * SN;
-  return (y_ref / sqrt(SN)).eval();
+  check_nonzero_size("unit_vector_constrain", "y", y);
+  return make_holder(
+      [&lp](const auto& y_ref) {
+        value_type_t<T1> SN = dot_self(y_ref);
+        check_positive_finite("unit_vector_constrain", "norm", SN);
+        lp -= 0.5 * SN;
+        return y_ref / sqrt(SN);
+      },
+      to_ref(y));
 }
 
 }  // namespace math
diff --git a/stan/math/prim/fun/unitspaced_array.hpp b/stan/math/prim/fun/unitspaced_array.hpp
index 809be55b528..4985ee2a5bb 100644
--- a/stan/math/prim/fun/unitspaced_array.hpp
+++ b/stan/math/prim/fun/unitspaced_array.hpp
@@ -24,10 +24,11 @@ inline std::vector<int> unitspaced_array(int low, int high) {
   static const char* function = "unitspaced_array";
   check_greater_or_equal(function, "high", high, low);
 
-  int K = fabs(high - low + 1);
-
-  Eigen::VectorXi v = Eigen::VectorXi::LinSpaced(K, low, high);
-  return {&v[0], &v[0] + K};
+  int K = std::abs(high - low + 1);
+  std::vector<int> result(K);
+  Eigen::Map<Eigen::VectorXi>(result.data(), K)
+      = Eigen::VectorXi::LinSpaced(K, low, high);
+  return result;
 }
 
 }  // namespace math
diff --git a/stan/math/prim/fun/zeros_row_vector.hpp b/stan/math/prim/fun/zeros_row_vector.hpp
index 1ba25e879fa..9b9f8df57cc 100644
--- a/stan/math/prim/fun/zeros_row_vector.hpp
+++ b/stan/math/prim/fun/zeros_row_vector.hpp
@@ -14,7 +14,7 @@ namespace math {
  * @return A row vector of size K with all elements initialised to 0.
  * @throw std::domain_error if K is negative.
  */
-inline Eigen::RowVectorXd zeros_row_vector(int K) {
+inline auto zeros_row_vector(int K) {
   check_nonnegative("zeros_row_vector", "size", K);
   return Eigen::RowVectorXd::Zero(K);
 }
diff --git a/stan/math/prim/fun/zeros_vector.hpp b/stan/math/prim/fun/zeros_vector.hpp
index 5e9acee53b5..024e97e93cf 100644
--- a/stan/math/prim/fun/zeros_vector.hpp
+++ b/stan/math/prim/fun/zeros_vector.hpp
@@ -14,7 +14,7 @@ namespace math {
  * @return A vector of size K with all elements initialised to 0.
  * @throw std::domain_error if K is negative.
  */
-inline Eigen::VectorXd zeros_vector(int K) {
+inline auto zeros_vector(int K) {
   check_nonnegative("zeros_vector", "size", K);
   return Eigen::VectorXd::Zero(K);
 }
diff --git a/stan/math/prim/functor/apply_scalar_binary.hpp b/stan/math/prim/functor/apply_scalar_binary.hpp
index c52a88586fd..f3a2c54ff15 100644
--- a/stan/math/prim/functor/apply_scalar_binary.hpp
+++ b/stan/math/prim/functor/apply_scalar_binary.hpp
@@ -55,7 +55,7 @@ template <typename T1, typename T2, typename F,
           require_all_eigen_t<T1, T2>* = nullptr>
 inline auto apply_scalar_binary(const T1& x, const T2& y, const F& f) {
   check_matching_dims("Binary function", "x", x, "y", y);
-  return x.binaryExpr(y, f).eval();
+  return x.binaryExpr(y, f);
 }
 
 /**
@@ -77,7 +77,7 @@ inline auto apply_scalar_binary(const T1& x, const T2& y, const F& f) {
   check_matching_sizes("Binary function", "x", x, "y", y);
   using int_vec_t = promote_scalar_t<value_type_t<T2>, plain_type_t<T1>>;
   Eigen::Map<const int_vec_t> y_map(y.data(), y.size());
-  return x.binaryExpr(y_map, f).eval();
+  return x.binaryExpr(y_map, f);
 }
 
 /**
@@ -99,7 +99,7 @@ inline auto apply_scalar_binary(const T1& x, const T2& y, const F& f) {
   check_matching_sizes("Binary function", "x", x, "y", y);
   using int_vec_t = promote_scalar_t<value_type_t<T1>, plain_type_t<T2>>;
   Eigen::Map<const int_vec_t> x_map(x.data(), x.size());
-  return x_map.binaryExpr(y, f).eval();
+  return x_map.binaryExpr(y, f);
 }
 
 /**
@@ -189,7 +189,7 @@ inline auto apply_scalar_binary(const T1& x, const T2& y, const F& f) {
 template <typename T1, typename T2, typename F, require_eigen_t<T1>* = nullptr,
           require_stan_scalar_t<T2>* = nullptr>
 inline auto apply_scalar_binary(const T1& x, const T2& y, const F& f) {
-  return x.unaryExpr([&f, &y](const auto& v) { return f(v, y); }).eval();
+  return x.unaryExpr([&f, &y](const auto& v) { return f(v, y); });
 }
 
 /**
@@ -213,7 +213,7 @@ inline auto apply_scalar_binary(const T1& x, const T2& y, const F& f) {
 template <typename T1, typename T2, typename F,
           require_stan_scalar_t<T1>* = nullptr, require_eigen_t<T2>* = nullptr>
 inline auto apply_scalar_binary(const T1& x, const T2& y, const F& f) {
-  return y.unaryExpr([&f, &x](const auto& v) { return f(x, v); }).eval();
+  return y.unaryExpr([&f, &x](const auto& v) { return f(x, v); });
 }
 
 /**
@@ -324,7 +324,7 @@ template <typename T1, typename T2, typename F,
           require_all_std_vector_vt<is_container, T1, T2>* = nullptr>
 inline auto apply_scalar_binary(const T1& x, const T2& y, const F& f) {
   check_matching_sizes("Binary function", "x", x, "y", y);
-  using T_return = decltype(apply_scalar_binary(x[0], y[0], f));
+  using T_return = plain_type_t<decltype(apply_scalar_binary(x[0], y[0], f))>;
   size_t y_size = y.size();
   std::vector<T_return> result(y_size);
   for (size_t i = 0; i < y_size; ++i) {
@@ -351,7 +351,7 @@ template <typename T1, typename T2, typename F,
           require_std_vector_vt<is_container, T1>* = nullptr,
           require_stan_scalar_t<T2>* = nullptr>
 inline auto apply_scalar_binary(const T1& x, const T2& y, const F& f) {
-  using T_return = decltype(apply_scalar_binary(x[0], y, f));
+  using T_return = plain_type_t<decltype(apply_scalar_binary(x[0], y, f))>;
   size_t x_size = x.size();
   std::vector<T_return> result(x_size);
   for (size_t i = 0; i < x_size; ++i) {
@@ -378,7 +378,7 @@ template <typename T1, typename T2, typename F,
           require_stan_scalar_t<T1>* = nullptr,
           require_std_vector_vt<is_container, T2>* = nullptr>
 inline auto apply_scalar_binary(const T1& x, const T2& y, const F& f) {
-  using T_return = decltype(apply_scalar_binary(x, y[0], f));
+  using T_return = plain_type_t<decltype(apply_scalar_binary(x, y[0], f))>;
   size_t y_size = y.size();
   std::vector<T_return> result(y_size);
   for (size_t i = 0; i < y_size; ++i) {
diff --git a/stan/math/prim/functor/apply_scalar_unary.hpp b/stan/math/prim/functor/apply_scalar_unary.hpp
index 53c5c651597..6a249dd92cd 100644
--- a/stan/math/prim/functor/apply_scalar_unary.hpp
+++ b/stan/math/prim/functor/apply_scalar_unary.hpp
@@ -5,6 +5,7 @@
 #include <stan/math/prim/meta/is_eigen.hpp>
 #include <stan/math/prim/meta/is_vector.hpp>
 #include <stan/math/prim/meta/is_vector_like.hpp>
+#include <stan/math/prim/meta/plain_type.hpp>
 #include <utility>
 #include <vector>
 
@@ -61,11 +62,8 @@ struct apply_scalar_unary<F, T, require_eigen_t<T>> {
    * by F to the specified matrix.
    */
   static inline auto apply(const T& x) {
-    return x
-        .unaryExpr([](scalar_t x) {
-          return apply_scalar_unary<F, scalar_t>::apply(x);
-        })
-        .eval();
+    return x.unaryExpr(
+        [](scalar_t x) { return apply_scalar_unary<F, scalar_t>::apply(x); });
   }
 
   /**
@@ -73,7 +71,7 @@ struct apply_scalar_unary<F, T, require_eigen_t<T>> {
    * expression template of type T.
    */
   using return_t = std::decay_t<decltype(
-      apply_scalar_unary<F, T>::apply(std::declval<T>()).eval())>;
+      apply_scalar_unary<F, T>::apply(std::declval<T>()))>;
 };
 
 /**
@@ -143,8 +141,8 @@ struct apply_scalar_unary<F, std::vector<T>> {
    * Return type, which is calculated recursively as a standard
    * vector of the return type of the contained type T.
    */
-  using return_t =
-      typename std::vector<typename apply_scalar_unary<F, T>::return_t>;
+  using return_t = typename std::vector<
+      plain_type_t<typename apply_scalar_unary<F, T>::return_t>>;
 
   /**
    * Apply the function specified by F elementwise to the
diff --git a/stan/math/prim/functor/apply_vector_unary.hpp b/stan/math/prim/functor/apply_vector_unary.hpp
index 022e97aaee6..2e2862e9081 100644
--- a/stan/math/prim/functor/apply_vector_unary.hpp
+++ b/stan/math/prim/functor/apply_vector_unary.hpp
@@ -12,7 +12,6 @@
 
 namespace stan {
 namespace math {
-
 // Forward declaration to allow specialisations
 template <typename T, typename Enable = void>
 struct apply_vector_unary {};
@@ -33,11 +32,7 @@ template <typename T>
 struct apply_vector_unary<T, require_eigen_t<T>> {
   /**
    * Member function for applying a functor to a vector and subsequently
-   * returning a vector. The returned objects are evaluated so that
-   * expression templates are not propagated.
-   *
-   * TODO(Andrew) Remove .eval() when rest of Stan library can take
-   * expressions as inputs.
+   * returning a vector.
    *
    * @tparam T Type of argument to which functor is applied.
    * @tparam F Type of functor to apply.
@@ -48,13 +43,38 @@ struct apply_vector_unary<T, require_eigen_t<T>> {
   template <typename F, typename T2 = T,
             require_t<is_eigen_matrix_base<plain_type_t<T2>>>* = nullptr>
   static inline auto apply(const T& x, const F& f) {
-    return f(x).matrix().eval();
+    return make_holder([](const auto& a) { return a.matrix().derived(); },
+                       f(x));
   }
 
   template <typename F, typename T2 = T,
             require_t<is_eigen_array<plain_type_t<T2>>>* = nullptr>
   static inline auto apply(const T& x, const F& f) {
-    return f(x).array().eval();
+    return make_holder([](const auto& a) { return a.array().derived(); }, f(x));
+  }
+
+  /**
+   * Member function for applying a functor to a vector and subsequently
+   * returning a vector. This is a variant of `apply` that does not construct
+   * `holder`, so it is up to the caller to ensure the returned expression is
+   * evaluated before `x` is destructed.
+   *
+   * @tparam T Type of argument to which functor is applied.
+   * @tparam F Type of functor to apply.
+   * @param x Eigen input to which operation is applied.
+   * @param f functor to apply to Eigen input.
+   * @return Eigen object with result of applying functor to input
+   */
+  template <typename F, typename T2 = T,
+            require_t<is_eigen_matrix_base<plain_type_t<T2>>>* = nullptr>
+  static inline auto apply_no_holder(const T& x, const F& f) {
+    return f(x).matrix().derived();
+  }
+
+  template <typename F, typename T2 = T,
+            require_t<is_eigen_array<plain_type_t<T2>>>* = nullptr>
+  static inline auto apply_no_holder(const T& x, const F& f) {
+    return f(x).array().derived();
   }
 
   /**
@@ -107,6 +127,22 @@ struct apply_vector_unary<T, require_std_vector_vt<is_stan_scalar, T>> {
     return result;
   }
 
+  /**
+   * Member function for applying a functor to each container in an std::vector
+   * and subsequently returning an std::vector of containers.
+   *
+   * @tparam T Type of argument to which functor is applied.
+   * @tparam F Type of functor to apply.
+   * @param x std::vector of containers to which operation is applied.
+   * @param f functor to apply to vector input.
+   * @return std::vector of containers with result of applying functor to
+   *         input.
+   */
+  template <typename F>
+  static inline auto apply_no_holder(const T& x, const F& f) {
+    return apply(x, f);
+  }
+
   /**
    * Member function for applying a functor to a vector and subsequently
    * returning a scalar.
@@ -151,13 +187,30 @@ struct apply_vector_unary<T, require_std_vector_vt<is_container, T>> {
   template <typename F>
   static inline auto apply(const T& x, const F& f) {
     size_t x_size = x.size();
-    using T_return = decltype(apply_vector_unary<T_vt>::apply(x[0], f));
+    using T_return
+        = plain_type_t<decltype(apply_vector_unary<T_vt>::apply(x[0], f))>;
     std::vector<T_return> result(x_size);
     for (size_t i = 0; i < x_size; ++i)
-      result[i] = apply_vector_unary<T_vt>::apply(x[i], f);
+      result[i] = apply_vector_unary<T_vt>::apply_no_holder(x[i], f);
     return result;
   }
 
+  /**
+   * Member function for applying a functor to each container in an std::vector
+   * and subsequently returning an std::vector of containers.
+   *
+   * @tparam T Type of argument to which functor is applied.
+   * @tparam F Type of functor to apply.
+   * @param x std::vector of containers to which operation is applied.
+   * @param f functor to apply to vector input.
+   * @return std::vector of containers with result of applying functor to
+   *         input.
+   */
+  template <typename F>
+  static inline auto apply_no_holder(const T& x, const F& f) {
+    return apply(x, f);
+  }
+
   /**
    * Member function for applying a functor to each container in an
    * std::vector and subsequently returning an std::vector of scalars.
diff --git a/stan/math/prim/meta/holder.hpp b/stan/math/prim/meta/holder.hpp
index 238cf854578..796caf8a541 100644
--- a/stan/math/prim/meta/holder.hpp
+++ b/stan/math/prim/meta/holder.hpp
@@ -2,6 +2,7 @@
 #define STAN_MATH_PRIM_FUN_HOLDER_HPP
 
 #include <stan/math/prim/meta/is_eigen.hpp>
+#include <stan/math/prim/meta/is_plain_type.hpp>
 #include <stan/math/prim/fun/Eigen.hpp>
 #include <memory>
 #include <type_traits>
@@ -248,9 +249,9 @@ T holder(T&& arg) {
 namespace internal {
 // the function holder_handle_element is also used in holder_cl
 /**
- * Handles single element (moving rvalues to heap) for construction of
- * `holder` or `holder_cl` from a functor. For lvalues just sets the `res`
- * pointer.
+ * Handles single element (moving rvalue non-expressions to heap) for
+ * construction of `holder` or `holder_cl` from a functor. For lvalues or
+ * expressions just sets the `res` pointer.
  * @tparam T type of the element
  * @param a element to handle
  * @param res resulting pointer to element
@@ -261,17 +262,34 @@ auto holder_handle_element(T& a, T*& res) {
   res = &a;
   return std::make_tuple();
 }
+template <typename T,
+          std::enable_if_t<!(Eigen::internal::traits<std::decay_t<T>>::Flags
+                             & Eigen::NestByRefBit)>* = nullptr>
+auto holder_handle_element(T&& a, std::remove_reference_t<T>*& res) {
+  res = &a;
+  return std::make_tuple();
+}
 
 /**
- * Handles single element (moving rvalues to heap) for construction of
- * `holder` or `holder_cl` from a functor. Rvalue is moved to heap and the
- * pointer to heap memory is assigned to res and returned in a tuple.
+ * Handles single element (moving rvalue non-expressions to heap) for
+ * construction of `holder` or `holder_cl` from a functor. Rvalue non-expression
+ * is moved to heap and the pointer to heap memory is assigned to res and
+ * returned in a tuple.
  * @tparam T type of the element
  * @param a element to handle
  * @param res resulting pointer to element
  * @return tuple of pointers allocated on heap (containing single pointer).
  */
-template <typename T, require_t<std::is_rvalue_reference<T&&>>* = nullptr>
+template <typename T, require_t<std::is_rvalue_reference<T&&>>* = nullptr,
+          std::enable_if_t<
+              static_cast<bool>(Eigen::internal::traits<std::decay_t<T>>::Flags&
+                                    Eigen::NestByRefBit)>* = nullptr>
+auto holder_handle_element(T&& a, T*& res) {
+  res = new T(std::move(a));
+  return std::make_tuple(res);
+}
+template <typename T, require_t<std::is_rvalue_reference<T&&>>* = nullptr,
+          require_not_eigen_t<T>* = nullptr>
 auto holder_handle_element(T&& a, T*& res) {
   res = new T(std::move(a));
   return std::make_tuple(res);
@@ -321,6 +339,7 @@ auto make_holder_impl(const F& func, std::index_sequence<Is...>,
  * This is similar to calling the functor with given arguments. Except that any
  * rvalue argument will be moved to heap first. The arguments moved to heap are
  * deleted once the expression is destructed.
+ *
  * @tparam F type of the functor
  * @tparam Args types of the arguments
  * @param func the functor
@@ -328,7 +347,7 @@ auto make_holder_impl(const F& func, std::index_sequence<Is...>,
  * @return `holder` referencing expression constructed by given functor
  */
 template <typename F, typename... Args,
-          require_eigen_t<
+          require_not_plain_type_t<
               decltype(std::declval<F>()(std::declval<Args&>()...))>* = nullptr>
 auto make_holder(const F& func, Args&&... args) {
   return internal::make_holder_impl(func,
@@ -336,6 +355,23 @@ auto make_holder(const F& func, Args&&... args) {
                                     std::forward<Args>(args)...);
 }
 
+/**
+ * Calls given function with given arguments. No `holder` is necessary if the
+ * function is not returning Eigen expression.
+ *
+ * @tparam F type of the functor
+ * @tparam Args types of the arguments
+ * @param func the functor
+ * @param args arguments for the functor
+ * @return `holder` referencing expression constructed by given functor
+ */
+template <typename F, typename... Args,
+          require_plain_type_t<
+              decltype(std::declval<F>()(std::declval<Args&>()...))>* = nullptr>
+auto make_holder(const F& func, Args&&... args) {
+  return func(std::forward<Args>(args)...);
+}
+
 }  // namespace math
 }  // namespace stan
 
diff --git a/stan/math/prim/meta/plain_type.hpp b/stan/math/prim/meta/plain_type.hpp
index ef742a5bdab..511ff730545 100644
--- a/stan/math/prim/meta/plain_type.hpp
+++ b/stan/math/prim/meta/plain_type.hpp
@@ -1,7 +1,6 @@
 #ifndef STAN_MATH_PRIM_META_PLAIN_TYPE_HPP
 #define STAN_MATH_PRIM_META_PLAIN_TYPE_HPP
 
-#include <stan/math/prim/meta/plain_type.hpp>
 #include <stan/math/prim/meta/is_eigen.hpp>
 #include <stan/math/prim/meta/is_detected.hpp>
 #include <type_traits>
diff --git a/stan/math/prim/prob/poisson_lpmf.hpp b/stan/math/prim/prob/poisson_lpmf.hpp
index 8fe7b769c94..f78420a5002 100644
--- a/stan/math/prim/prob/poisson_lpmf.hpp
+++ b/stan/math/prim/prob/poisson_lpmf.hpp
@@ -11,6 +11,7 @@
 #include <stan/math/prim/fun/promote_scalar.hpp>
 #include <stan/math/prim/fun/size.hpp>
 #include <stan/math/prim/fun/size_zero.hpp>
+#include <stan/math/prim/fun/sum.hpp>
 #include <stan/math/prim/fun/value_of.hpp>
 #include <stan/math/prim/functor/operands_and_partials.hpp>
 
@@ -66,7 +67,7 @@ return_type_t<T_rate> poisson_lpmf(const T_n& n, const T_rate& lambda) {
 
   operands_and_partials<T_lambda_ref> ops_partials(lambda_ref);
 
-  T_partials_return logp = sum(multiply_log(n_val, lambda_val));
+  T_partials_return logp = stan::math::sum(multiply_log(n_val, lambda_val));
   if (include_summand<propto, T_rate>::value) {
     logp -= sum(lambda_val) * N / size(lambda);
   }
diff --git a/stan/math/rev.hpp b/stan/math/rev.hpp
index 3b56ffe3abc..63ef295d138 100644
--- a/stan/math/rev.hpp
+++ b/stan/math/rev.hpp
@@ -5,10 +5,9 @@
 
 #include <stan/math/rev/core.hpp>
 #include <stan/math/rev/meta.hpp>
-
-#include <stan/math/prim.hpp>
-
 #include <stan/math/rev/fun.hpp>
 #include <stan/math/rev/functor.hpp>
 
+#include <stan/math/prim.hpp>
+
 #endif
diff --git a/stan/math/rev/fun.hpp b/stan/math/rev/fun.hpp
index d53279c4588..4b9de3df0ba 100644
--- a/stan/math/rev/fun.hpp
+++ b/stan/math/rev/fun.hpp
@@ -3,7 +3,10 @@
 
 #include <stan/math/prim/fun/Eigen.hpp>
 
-#include <stan/math/prim.hpp>
+#include <stan/math/prim/err.hpp>
+#include <stan/math/prim/fun.hpp>
+#include <stan/math/prim/functor.hpp>
+#include <stan/math/prim/meta.hpp>
 
 #include <stan/math/rev/fun/LDLT_alloc.hpp>
 #include <stan/math/rev/fun/LDLT_factor.hpp>
diff --git a/stan/math/rev/fun/multiply.hpp b/stan/math/rev/fun/multiply.hpp
index 7ded4e0b147..4f7ebee9cd9 100644
--- a/stan/math/rev/fun/multiply.hpp
+++ b/stan/math/rev/fun/multiply.hpp
@@ -4,7 +4,7 @@
 #include <stan/math/rev/meta.hpp>
 #include <stan/math/rev/core.hpp>
 #include <stan/math/rev/core/typedefs.hpp>
-#include <stan/math/prim.hpp>
+#include <stan/math/prim/fun.hpp>
 #include <type_traits>
 
 namespace stan {
diff --git a/stan/math/rev/functor/apply_vector_unary.hpp b/stan/math/rev/functor/apply_vector_unary.hpp
index d244047108b..612b767f969 100644
--- a/stan/math/rev/functor/apply_vector_unary.hpp
+++ b/stan/math/rev/functor/apply_vector_unary.hpp
@@ -29,6 +29,21 @@ struct apply_vector_unary<T, require_var_matrix_t<T>> {
     return f(x);
   }
 
+  /**
+   * Member function for applying a functor to a `var_value<T>` and
+   * subsequently returning a `var_value<T>`.
+   *
+   * @tparam T Type of argument to which functor is applied.
+   * @tparam F Type of functor to apply.
+   * @param x input to which operation is applied.
+   * @param f functor to apply to Eigen input.
+   * @return object with result of applying functor to input
+   */
+  template <typename F>
+  static inline T apply_no_holder(const T& x, const F& f) {
+    return f(x);
+  }
+
   /**
    * Member function for applying a functor to a `var_value<T>` and
    * subsequently returning a var. The reduction to a var needs
diff --git a/test/unit/math/prim/fun/binary_scalar_tester.hpp b/test/unit/math/prim/fun/binary_scalar_tester.hpp
index 51fbb930b3b..747b948f99f 100644
--- a/test/unit/math/prim/fun/binary_scalar_tester.hpp
+++ b/test/unit/math/prim/fun/binary_scalar_tester.hpp
@@ -21,9 +21,9 @@ namespace test {
 template <typename F, typename T1, typename T2,
           require_all_not_vector_t<T1, T2>* = nullptr>
 void binary_scalar_tester_impl(const F& f, const T1& x, const T2& y) {
-  auto vec_vec = f(x, y);
-  auto vec_scal = f(x, y(0));
-  auto scal_vec = f(x(0), y);
+  auto vec_vec = math::eval(f(x, y));
+  auto vec_scal = math::eval(f(x, y(0)));
+  auto scal_vec = math::eval(f(x(0), y));
   for (int i = 0; i < x.size(); ++i) {
     EXPECT_FLOAT_EQ(f(x(i), y(i)), vec_vec(i));
     EXPECT_FLOAT_EQ(f(x(i), y(0)), vec_scal(i));
@@ -91,9 +91,9 @@ void binary_scalar_tester_impl(const F& f, const T1& x, const T2& y) {
 template <typename F, typename T1, typename T2,
           require_all_vector_t<T1, T2>* = nullptr>
 void binary_scalar_tester_impl(const F& f, const T1& x, const T2& y) {
-  auto vec_vec = f(x, y);
-  auto vec_scal = f(x, y[0]);
-  auto scal_vec = f(x[0], y);
+  auto vec_vec = math::eval(f(x, y));
+  auto vec_scal = math::eval(f(x, y[0]));
+  auto scal_vec = math::eval(f(x[0], y));
   for (int i = 0; i < x.size(); ++i) {
     EXPECT_FLOAT_EQ(f(x[i], y[i]), vec_vec[i]);
     EXPECT_FLOAT_EQ(f(x[i], y[0]), vec_scal[i]);
@@ -163,7 +163,7 @@ template <typename F, typename T1, typename T2,
           require_eigen_matrix_dynamic_t<T1>* = nullptr,
           require_std_vector_t<T2>* = nullptr>
 void binary_scalar_tester_impl(const F& f, const T1& x, const T2& y) {
-  auto vec_vec = f(x, y);
+  auto vec_vec = math::eval(f(x, y));
   for (int r = 0; r < x.rows(); ++r) {
     for (int c = 0; c < x.cols(); ++c) {
       EXPECT_FLOAT_EQ(f(x(r, c), y[r][c]), vec_vec(r, c));
@@ -216,7 +216,7 @@ template <typename F, typename T1, typename T2,
           require_std_vector_t<T1>* = nullptr,
           require_eigen_matrix_dynamic_t<T2>* = nullptr>
 void binary_scalar_tester_impl(const F& f, const T1& x, const T2& y) {
-  auto vec_vec = f(x, y);
+  auto vec_vec = math::eval(f(x, y));
   for (int r = 0; r < y.rows(); ++r) {
     for (int c = 0; c < y.cols(); ++c) {
       EXPECT_FLOAT_EQ(f(x[r][c], y(r, c)), vec_vec(r, c));
diff --git a/test/unit/math/serializer.hpp b/test/unit/math/serializer.hpp
index 38b46148847..81258076cba 100644
--- a/test/unit/math/serializer.hpp
+++ b/test/unit/math/serializer.hpp
@@ -194,7 +194,7 @@ struct serializer {
    * @tparam U type of specified scalar; must be assignable to T
    * @param x scalar to serialize
    */
-  template <typename U>
+  template <typename U, require_stan_scalar_t<U>* = nullptr>
   void write(const U& x) {
     vals_.push_back(x);
   }
@@ -231,10 +231,11 @@ struct serializer {
    * @tparam C column specification of Eigen container
    * @param x Eigen container to serialize.
    */
-  template <typename U, int R, int C>
-  void write(const Eigen::Matrix<U, R, C>& x) {
-    for (int i = 0; i < x.size(); ++i)
-      write(x(i));
+  template <typename U, require_eigen_t<U>* = nullptr>
+  void write(const U& x) {
+    for (int i = 0; i < x.size(); ++i) {
+      write(x.coeff(i));
+    }
   }
 
   /**