diff --git a/stan/math/prim/fun.hpp b/stan/math/prim/fun.hpp
index a09499863d4..8ef4204c81b 100644
--- a/stan/math/prim/fun.hpp
+++ b/stan/math/prim/fun.hpp
@@ -353,6 +353,6 @@
 #include <stan/math/prim/fun/zeros_int_array.hpp>
 #include <stan/math/prim/fun/zeros_row_vector.hpp>
 #include <stan/math/prim/fun/zeros_vector.hpp>
-#include <stan/math/prim/fun/zip.hpp>
+#include <stan/math/prim/fun/zip_index.hpp>
 
 #endif
diff --git a/stan/math/prim/fun/zip.hpp b/stan/math/prim/fun/zip_index.hpp
similarity index 64%
rename from stan/math/prim/fun/zip.hpp
rename to stan/math/prim/fun/zip_index.hpp
index 4c1903cd2bc..e9c93fdf0cd 100644
--- a/stan/math/prim/fun/zip.hpp
+++ b/stan/math/prim/fun/zip_index.hpp
@@ -1,5 +1,5 @@
-#ifndef STAN_MATH_PRIM_FUN_ZIP_HPP
-#define STAN_MATH_PRIM_FUN_ZIP_HPP
+#ifndef STAN_MATH_PRIM_FUN_ZIP_INDEX_HPP
+#define STAN_MATH_PRIM_FUN_ZIP_INDEX_HPP
 
 #include <stan/math/prim/meta.hpp>
 #include <stan/math/prim/err.hpp>
@@ -10,9 +10,8 @@ namespace stan {
 namespace math {
 
 /**
- * Return a vector of matrix values obtained by zipping
- * two N-dimensional integer vectors to form tuples of
- * corresponding elements.
+ * Return the column vector of elements of x selected
+ * by paired (zipped) 1-based row and column indices
  *
  * @tparam EigMat type of the matrix
  * @param x input matrix
@@ -25,9 +24,9 @@ namespace math {
 template <typename EigMat, typename IdxRows, typename IdxCols,
           require_eigen_matrix_dynamic_t<EigMat>* = nullptr,
           require_all_vector_t<IdxRows, IdxCols>* = nullptr>
-inline auto zip(EigMat&& x, IdxRows&& idx_row, IdxCols&& idx_col) {
-  check_size_match("zip", "size of idx_row", idx_row.size(), "size of idx_col",
-                   idx_col.size());
+inline auto zip_index(EigMat&& x, IdxRows&& idx_row, IdxCols&& idx_col) {
+  check_size_match("zip_index", "size of idx_row", idx_row.size(),
+                   "size of idx_col", idx_col.size());
   return make_holder(
       [](auto&& x_, auto&& idx_row_, auto&& idx_col_) {
         using map_t = Eigen::Map<const Eigen::Array<int, Eigen::Dynamic, 1>>;
@@ -35,10 +34,14 @@ inline auto zip(EigMat&& x, IdxRows&& idx_row, IdxCols&& idx_col) {
         const map_t cols(idx_col_.data(), idx_col_.size());
     // If the user turns of range checks do not pay for min and max sweeps
 #ifndef STAN_NO_RANGE_CHECKS
-        check_range("zip", "minimum row index", x_.rows(), rows.minCoeff());
-        check_range("zip", "maximum row index", x_.rows(), rows.maxCoeff());
-        check_range("zip", "minimum column index", x_.cols(), cols.minCoeff());
-        check_range("zip", "maximum column index", x_.cols(), cols.maxCoeff());
+        check_range("zip_index", "minimum row index", x_.rows(),
+                    rows.minCoeff());
+        check_range("zip_index", "maximum row index", x_.rows(),
+                    rows.maxCoeff());
+        check_range("zip_index", "minimum column index", x_.cols(),
+                    cols.minCoeff());
+        check_range("zip_index", "maximum column index", x_.cols(),
+                    cols.maxCoeff());
 #endif
         const auto linear_idx = (rows.cast<Eigen::Index>() - 1)
                                 + (cols.cast<Eigen::Index>() - 1) * x_.rows();
diff --git a/stan/math/rev/fun.hpp b/stan/math/rev/fun.hpp
index d6f02c3f13f..e5ef938a302 100644
--- a/stan/math/rev/fun.hpp
+++ b/stan/math/rev/fun.hpp
@@ -183,7 +183,7 @@
 #include <stan/math/rev/fun/value_of.hpp>
 #include <stan/math/rev/fun/value_of_rec.hpp>
 #include <stan/math/rev/fun/variance.hpp>
-#include <stan/math/rev/fun/zip.hpp>
+#include <stan/math/rev/fun/zip_index.hpp>
 
 #include <stan/math/prim/constraint.hpp>
 #include <stan/math/prim/err.hpp>
diff --git a/stan/math/rev/fun/zip.hpp b/stan/math/rev/fun/zip_index.hpp
similarity index 71%
rename from stan/math/rev/fun/zip.hpp
rename to stan/math/rev/fun/zip_index.hpp
index 9214cb6ce25..0c137e7333e 100644
--- a/stan/math/rev/fun/zip.hpp
+++ b/stan/math/rev/fun/zip_index.hpp
@@ -1,11 +1,11 @@
-#ifndef STAN_MATH_REV_FUN_ZIP_HPP
-#define STAN_MATH_REV_FUN_ZIP_HPP
+#ifndef STAN_MATH_REV_FUN_ZIP_INDEX_HPP
+#define STAN_MATH_REV_FUN_ZIP_INDEX_HPP
 
 #include <stan/math/rev/meta.hpp>
 #include <stan/math/rev/core.hpp>
 #include <stan/math/rev/fun/to_arena.hpp>
 #include <stan/math/prim/err/check_size_match.hpp>
-#include <stan/math/prim/fun/zip.hpp>
+#include <stan/math/prim/fun/zip_index.hpp>
 #include <vector>
 
 namespace stan::math {
@@ -29,18 +29,19 @@ namespace stan::math {
 template <typename EigMat, typename IdxRows, typename IdxCols,
           require_eigen_matrix_dynamic_t<EigMat>* = nullptr,
           require_all_std_vector_t<IdxRows, IdxCols>* = nullptr>
-inline var_value<Eigen::VectorXd> zip(const var_value<EigMat>& x,
-                                      IdxRows&& idx_row, IdxCols&& idx_col) {
-  check_size_match("zip", "size of idx_row", idx_row.size(), "size of idx_col",
-                   idx_col.size());
+inline var_value<Eigen::VectorXd> zip_index(const var_value<EigMat>& x,
+                                            IdxRows&& idx_row,
+                                            IdxCols&& idx_col) {
+  check_size_match("zip_index", "size of idx_row", idx_row.size(),
+                   "size of idx_col", idx_col.size());
   if (idx_row.empty()) {
     return var_value<Eigen::VectorXd>(Eigen::VectorXd(0));
   }
   auto rows = to_arena(std::forward<IdxRows>(idx_row));
   auto cols = to_arena(std::forward<IdxCols>(idx_col));
-  return make_callback_var(zip(x.val(), rows, cols),
+  return make_callback_var(zip_index(x.val(), rows, cols),
                            [x, rows, cols](auto& vi) mutable {
-                             zip(x.adj(), rows, cols) += vi.adj();
+                             zip_index(x.adj(), rows, cols) += vi.adj();
                            });
 }
 
diff --git a/test/unit/math/prim/fun/zip_test.cpp b/test/unit/math/prim/fun/zip_index_test.cpp
similarity index 69%
rename from test/unit/math/prim/fun/zip_test.cpp
rename to test/unit/math/prim/fun/zip_index_test.cpp
index 79c489348e6..31b51f5304a 100644
--- a/test/unit/math/prim/fun/zip_test.cpp
+++ b/test/unit/math/prim/fun/zip_index_test.cpp
@@ -3,7 +3,7 @@
 #include <gtest/gtest.h>
 #include <vector>
 
-TEST(MathFunctions, zip) {
+TEST(MathFunctions, zip_index) {
   const int J = 5;  // Number of locations
   const int K = 4;  // Number of time categories
 
@@ -15,7 +15,7 @@ TEST(MathFunctions, zip) {
 
   std::vector<int> loc = {1, 1, 5, 2, 3, 3, 4, 4, 5};
 
-  Eigen::VectorXd zipped = stan::math::zip(delta, loc, time);
+  Eigen::VectorXd zipped = stan::math::zip_index(delta, loc, time);
 
   EXPECT_EQ(zipped.size(), loc.size());
   EXPECT_DOUBLE_EQ(zipped[0], delta(0, 0));
@@ -29,7 +29,7 @@ TEST(MathFunctions, zip) {
   EXPECT_DOUBLE_EQ(zipped[8], delta(4, 0));
 }
 
-TEST(MathFunctions, zip_errors) {
+TEST(MathFunctions, zip_index_errors) {
   Eigen::MatrixXd delta(5, 4);  // 5 rows, 4 columns
   delta.setZero();
 
@@ -37,25 +37,27 @@ TEST(MathFunctions, zip_errors) {
   std::vector<int> cols = {1, 2, 3};
 
   // valid input does not throw
-  EXPECT_NO_THROW(stan::math::zip(delta, rows, cols));
+  EXPECT_NO_THROW(stan::math::zip_index(delta, rows, cols));
 
   // index vectors of different sizes
   std::vector<int> cols_short = {1, 2};
-  EXPECT_THROW(stan::math::zip(delta, rows, cols_short), std::invalid_argument);
+  EXPECT_THROW(stan::math::zip_index(delta, rows, cols_short),
+               std::invalid_argument);
 
   // row index too large (matrix has 5 rows)
   std::vector<int> rows_big = {1, 6, 3};
-  EXPECT_THROW(stan::math::zip(delta, rows_big, cols), std::out_of_range);
+  EXPECT_THROW(stan::math::zip_index(delta, rows_big, cols), std::out_of_range);
 
   // row index too small (indices are 1-based)
   std::vector<int> rows_zero = {0, 2, 3};
-  EXPECT_THROW(stan::math::zip(delta, rows_zero, cols), std::out_of_range);
+  EXPECT_THROW(stan::math::zip_index(delta, rows_zero, cols),
+               std::out_of_range);
 
   // column index too large (matrix has 4 columns)
   std::vector<int> cols_big = {1, 5, 3};
-  EXPECT_THROW(stan::math::zip(delta, rows, cols_big), std::out_of_range);
+  EXPECT_THROW(stan::math::zip_index(delta, rows, cols_big), std::out_of_range);
 
   // column index too small
   std::vector<int> cols_neg = {1, -1, 3};
-  EXPECT_THROW(stan::math::zip(delta, rows, cols_neg), std::out_of_range);
+  EXPECT_THROW(stan::math::zip_index(delta, rows, cols_neg), std::out_of_range);
 }
diff --git a/test/unit/math/rev/fun/zip_test.cpp b/test/unit/math/rev/fun/zip_index_test.cpp
similarity index 74%
rename from test/unit/math/rev/fun/zip_test.cpp
rename to test/unit/math/rev/fun/zip_index_test.cpp
index f0e28e417c7..24316ca300e 100644
--- a/test/unit/math/rev/fun/zip_test.cpp
+++ b/test/unit/math/rev/fun/zip_index_test.cpp
@@ -3,13 +3,13 @@
 #include <gtest/gtest.h>
 #include <vector>
 
-TEST_F(AgradRev, zip_values_and_gradients) {
+TEST_F(AgradRev, zip_index_values_and_gradients) {
   Eigen::MatrixXd values(3, 2);
   values << 11, 12, 21, 22, 31, 32;
   stan::math::var_value<Eigen::MatrixXd> x = values;
   const std::vector<int> rows{3, 1, 3, 2, 1};
   const std::vector<int> cols{2, 1, 2, 1, 2};
-  const auto y = stan::math::zip(x, rows, cols);
+  const auto y = stan::math::zip_index(x, rows, cols);
 
   EXPECT_EQ(y.rows(), 5);
   EXPECT_EQ(y.cols(), 1);
@@ -36,12 +36,12 @@ TEST_F(AgradRev, zip_index_lifetime) {
   auto y = [&x]() {
     std::vector<int> rows{1, 2, 1};
     std::vector<int> cols{3, 1, 3};
-    auto result = stan::math::zip(x, std::move(rows), std::move(cols));
+    auto result = stan::math::zip_index(x, std::move(rows), std::move(cols));
     rows.assign(3, 2);
     cols.assign(3, 2);
     return result;
   }();
-  auto z = stan::math::zip(x, std::vector<int>{2}, std::vector<int>{2});
+  auto z = stan::math::zip_index(x, std::vector<int>{2}, std::vector<int>{2});
   stan::math::var objective = stan::math::sum(y) + stan::math::sum(z);
   objective.grad();
   EXPECT_DOUBLE_EQ(x.adj()(0, 0), 0);
@@ -52,12 +52,12 @@ TEST_F(AgradRev, zip_index_lifetime) {
   EXPECT_DOUBLE_EQ(x.adj()(1, 2), 0);
 }
 
-TEST_F(AgradRev, zip_matrix_view) {
+TEST_F(AgradRev, zip_index_matrix_view) {
   Eigen::MatrixXd values(3, 3);
   values << 11, 12, 13, 21, 22, 23, 31, 32, 33;
   stan::math::var_value<Eigen::MatrixXd> x = values;
-  auto y = stan::math::zip(x.block(1, 1, 2, 2), std::vector<int>{2, 1, 2},
-                           std::vector<int>{1, 2, 1});
+  auto y = stan::math::zip_index(x.block(1, 1, 2, 2), std::vector<int>{2, 1, 2},
+                                 std::vector<int>{1, 2, 1});
   EXPECT_DOUBLE_EQ(y.val()(0), 32);
   EXPECT_DOUBLE_EQ(y.val()(1), 23);
   EXPECT_DOUBLE_EQ(y.val()(2), 32);
@@ -70,14 +70,14 @@ TEST_F(AgradRev, zip_matrix_view) {
   }
 }
 
-TEST_F(AgradRev, zip_row_major) {
+TEST_F(AgradRev, zip_index_row_major) {
   using matrix_t
       = Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>;
   matrix_t values(2, 3);
   values << 11, 12, 13, 21, 22, 23;
   stan::math::var_value<matrix_t> x = values;
-  auto y = stan::math::zip(x, std::vector<int>{2, 1, 2},
-                           std::vector<int>{1, 3, 1});
+  auto y = stan::math::zip_index(x, std::vector<int>{2, 1, 2},
+                                 std::vector<int>{1, 3, 1});
   EXPECT_DOUBLE_EQ(y.val()(0), 21);
   EXPECT_DOUBLE_EQ(y.val()(1), 13);
   EXPECT_DOUBLE_EQ(y.val()(2), 21);
@@ -90,13 +90,13 @@ TEST_F(AgradRev, zip_row_major) {
   }
 }
 
-TEST_F(AgradRev, zip_empty) {
+TEST_F(AgradRev, zip_index_empty) {
   const std::vector<int> indices;
   for (int rows : {0, 2}) {
     for (int cols : {0, 3}) {
       stan::math::var_value<Eigen::MatrixXd> x
           = Eigen::MatrixXd::Ones(rows, cols);
-      auto y = stan::math::zip(x, indices, indices);
+      auto y = stan::math::zip_index(x, indices, indices);
       EXPECT_EQ(y.rows(), 0);
       EXPECT_EQ(y.cols(), 1);
       stan::math::sum(y).grad();
@@ -105,18 +105,18 @@ TEST_F(AgradRev, zip_empty) {
   }
 }
 
-TEST_F(AgradRev, zip_errors) {
+TEST_F(AgradRev, zip_index_errors) {
   stan::math::var_value<Eigen::MatrixXd> x = Eigen::MatrixXd::Zero(3, 2);
   const std::vector<int> valid{1, 2};
   const std::vector<int> empty;
-  EXPECT_THROW(stan::math::zip(x, valid, empty), std::invalid_argument);
-  EXPECT_THROW(stan::math::zip(x, empty, valid), std::invalid_argument);
+  EXPECT_THROW(stan::math::zip_index(x, valid, empty), std::invalid_argument);
+  EXPECT_THROW(stan::math::zip_index(x, empty, valid), std::invalid_argument);
   for (int row : {-1, 0, 4}) {
-    EXPECT_THROW(stan::math::zip(x, std::vector<int>{1, row}, valid),
+    EXPECT_THROW(stan::math::zip_index(x, std::vector<int>{1, row}, valid),
                  std::out_of_range);
   }
   for (int col : {-1, 0, 3}) {
-    EXPECT_THROW(stan::math::zip(x, valid, std::vector<int>{1, col}),
+    EXPECT_THROW(stan::math::zip_index(x, valid, std::vector<int>{1, col}),
                  std::out_of_range);
   }
   for (int rows : {0, 2}) {
@@ -126,7 +126,8 @@ TEST_F(AgradRev, zip_errors) {
       }
       stan::math::var_value<Eigen::MatrixXd> empty_x
           = Eigen::MatrixXd::Zero(rows, cols);
-      EXPECT_THROW(stan::math::zip(empty_x, valid, valid), std::out_of_range);
+      EXPECT_THROW(stan::math::zip_index(empty_x, valid, valid),
+                   std::out_of_range);
     }
   }
 }