diff --git a/src/stan/services/optimize/bfgs.hpp b/src/stan/services/optimize/bfgs.hpp
index 37cd0a58973..6f74388eb11 100644
--- a/src/stan/services/optimize/bfgs.hpp
+++ b/src/stan/services/optimize/bfgs.hpp
@@ -83,6 +83,7 @@ int bfgs(Model& model, const stan::io::var_context& init,
   bfgs._conv_opts.maxIts = num_iterations;
 
   double lp = bfgs.logp();
+  int ret = 0;
 
   std::stringstream initial_msg;
   initial_msg << "Initial log joint probability = " << lp;
@@ -90,6 +91,7 @@ int bfgs(Model& model, const stan::io::var_context& init,
 
   std::vector<std::string> names;
   names.push_back("lp__");
+  names.push_back("converged__");
   model.constrained_param_names(names, true, true);
   parameter_writer(names);
 
@@ -109,10 +111,9 @@ int bfgs(Model& model, const stan::io::var_context& init,
     if (msg.str().length() > 0)
       logger.info(msg);
 
-    values.insert(values.begin(), lp);
+    values.insert(values.begin(), {lp, static_cast<double>(ret)});
     parameter_writer(values);
   }
-  int ret = 0;
 
   try {
     while (ret == 0) {
@@ -168,7 +169,7 @@ int bfgs(Model& model, const stan::io::var_context& init,
         if (msg.str().length() > 0)
           logger.info(msg);
 
-        values.insert(values.begin(), lp);
+        values.insert(values.begin(), {lp, static_cast<double>(ret)});
         parameter_writer(values);
       }
     }
@@ -192,7 +193,7 @@ int bfgs(Model& model, const stan::io::var_context& init,
     }
     if (msg.str().length() > 0)
       logger.info(msg);
-    values.insert(values.begin(), lp);
+    values.insert(values.begin(), {lp, static_cast<double>(ret)});
     parameter_writer(values);
   }
 
@@ -201,10 +202,12 @@ int bfgs(Model& model, const stan::io::var_context& init,
   if (ret >= 0) {
     logger.info("Optimization terminated normally: ");
     logger.info("  " + error_string);
+    parameter_writer("Optimization terminated normally: " + error_string);
     return_code = error_codes::OK;
   } else {
     logger.error("Optimization terminated with error: ");
     logger.error("  " + error_string);
+    parameter_writer("Optimization terminated with error: " + error_string);
     return_code = error_codes::SOFTWARE;
   }
 
diff --git a/src/stan/services/optimize/lbfgs.hpp b/src/stan/services/optimize/lbfgs.hpp
index bf50e9088f6..b38f79c1c23 100644
--- a/src/stan/services/optimize/lbfgs.hpp
+++ b/src/stan/services/optimize/lbfgs.hpp
@@ -87,6 +87,7 @@ int lbfgs(Model& model, const stan::io::var_context& init,
   lbfgs._conv_opts.maxIts = num_iterations;
 
   double lp = lbfgs.logp();
+  int ret = 0;
 
   std::stringstream initial_msg;
   initial_msg << "Initial log joint probability = " << lp;
@@ -94,6 +95,7 @@ int lbfgs(Model& model, const stan::io::var_context& init,
 
   std::vector<std::string> names;
   names.push_back("lp__");
+  names.push_back("converged__");
   model.constrained_param_names(names, true, true);
   parameter_writer(names);
 
@@ -104,10 +106,9 @@ int lbfgs(Model& model, const stan::io::var_context& init,
     if (msg.str().length() > 0)
       logger.info(msg);
 
-    values.insert(values.begin(), lp);
+    values.insert(values.begin(), {lp, static_cast<double>(ret)});
     parameter_writer(values);
   }
-  int ret = 0;
 
   try {
     while (ret == 0) {
@@ -161,7 +162,7 @@ int lbfgs(Model& model, const stan::io::var_context& init,
         if (msg.str().length() > 0)
           logger.info(msg);
 
-        values.insert(values.begin(), lp);
+        values.insert(values.begin(), {lp, static_cast<double>(ret)});
         parameter_writer(values);
       }
     }
@@ -186,7 +187,7 @@ int lbfgs(Model& model, const stan::io::var_context& init,
     if (msg.str().length() > 0)
       logger.info(msg);
 
-    values.insert(values.begin(), lp);
+    values.insert(values.begin(), {lp, static_cast<double>(ret)});
     parameter_writer(values);
   }
 
@@ -196,10 +197,12 @@ int lbfgs(Model& model, const stan::io::var_context& init,
   if (ret >= 0) {
     logger.info("Optimization terminated normally: ");
     logger.info("  " + error_string);
+    parameter_writer("Optimization terminated normally: " + error_string);
     return_code = error_codes::OK;
   } else {
     logger.error("Optimization terminated with error: ");
     logger.error("  " + error_string);
+    parameter_writer("Optimization terminated with error: " + error_string);
     return_code = error_codes::SOFTWARE;
   }
 
diff --git a/src/stan/services/optimize/newton.hpp b/src/stan/services/optimize/newton.hpp
index db64f6e46c3..b74aa282e9d 100644
--- a/src/stan/services/optimize/newton.hpp
+++ b/src/stan/services/optimize/newton.hpp
@@ -85,7 +85,9 @@ int newton(Model& model, const stan::io::var_context& init,
   logger.info(msg);
 
   std::vector<std::string> names;
+
   names.push_back("lp__");
+  names.push_back("converged__");
   model.constrained_param_names(names, true, true);
   parameter_writer(names);
 
@@ -97,7 +99,7 @@ int newton(Model& model, const stan::io::var_context& init,
       model.write_array(rng, cont_vector, disc_vector, values, true, true, &ss);
       if (ss.str().length() > 0)
         logger.info(ss);
-      values.insert(values.begin(), lp);
+      values.insert(values.begin(), {lp, 0});
       parameter_writer(values);
     }
     interrupt();
@@ -121,7 +123,7 @@ int newton(Model& model, const stan::io::var_context& init,
     model.write_array(rng, cont_vector, disc_vector, values, true, true, &ss);
     if (ss.str().length() > 0)
       logger.info(ss);
-    values.insert(values.begin(), lp);
+    values.insert(values.begin(), {lp, 0});
     parameter_writer(values);
   }
   return error_codes::OK;
diff --git a/src/test/unit/services/optimize/bfgs_jacobian_test.cpp b/src/test/unit/services/optimize/bfgs_jacobian_test.cpp
index ea71306a7f6..a92ec7e7dab 100644
--- a/src/test/unit/services/optimize/bfgs_jacobian_test.cpp
+++ b/src/test/unit/services/optimize/bfgs_jacobian_test.cpp
@@ -34,10 +34,11 @@ TEST_F(ServicesOptimize, withJacobian) {
   EXPECT_TRUE(logger.find("Optimization terminated normally: "));
   EXPECT_FLOAT_EQ(return_code, 0);
 
-  ASSERT_EQ(2, parameter.names_.size());
+  ASSERT_EQ(3, parameter.names_.size());
   EXPECT_EQ("lp__", parameter.names_[0]);
-  EXPECT_EQ("sigma", parameter.names_[1]);
-  EXPECT_NEAR((3 + std::sqrt(13)) / 2, parameter.states_.back()[1], 0.0001);
+  EXPECT_EQ("converged__", parameter.names_[1]);
+  EXPECT_EQ("sigma", parameter.names_[2]);
+  EXPECT_NEAR((3 + std::sqrt(13)) / 2, parameter.states_.back()[2], 0.0001);
   EXPECT_GT(interrupt.call_count(), 0);
 }
 
@@ -58,9 +59,10 @@ TEST_F(ServicesOptimize, withoutJacobian) {
   EXPECT_TRUE(logger.find("Optimization terminated normally: "));
   EXPECT_FLOAT_EQ(return_code, 0);
 
-  ASSERT_EQ(2, parameter.names_.size());
+  ASSERT_EQ(3, parameter.names_.size());
   EXPECT_EQ("lp__", parameter.names_[0]);
-  EXPECT_EQ("sigma", parameter.names_[1]);
-  EXPECT_NEAR(3, parameter.states_.back()[1], 0.0001);
+  EXPECT_EQ("converged__", parameter.names_[1]);
+  EXPECT_EQ("sigma", parameter.names_[2]);
+  EXPECT_NEAR(3, parameter.states_.back()[2], 0.0001);
   EXPECT_GT(interrupt.call_count(), 1);
 }
diff --git a/src/test/unit/services/optimize/bfgs_test.cpp b/src/test/unit/services/optimize/bfgs_test.cpp
index 6f74ffaae15..ac376fe1910 100644
--- a/src/test/unit/services/optimize/bfgs_test.cpp
+++ b/src/test/unit/services/optimize/bfgs_test.cpp
@@ -40,20 +40,21 @@ TEST_F(ServicesOptimize, rosenbrock) {
 
   EXPECT_EQ("0,0\n", init_ss.str());
 
-  ASSERT_EQ(3, parameter.names_.size());
+  ASSERT_EQ(4, parameter.names_.size());
   EXPECT_EQ("lp__", parameter.names_[0]);
-  EXPECT_EQ("x", parameter.names_[1]);
-  EXPECT_EQ("y", parameter.names_[2]);
+  EXPECT_EQ("converged__", parameter.names_[1]);
+  EXPECT_EQ("x", parameter.names_[2]);
+  EXPECT_EQ("y", parameter.names_[3]);
 
   EXPECT_EQ(20, parameter.states_.size());
-  EXPECT_FLOAT_EQ(0, parameter.states_.front()[1])
-      << "initial value should be (0, 0)";
   EXPECT_FLOAT_EQ(0, parameter.states_.front()[2])
       << "initial value should be (0, 0)";
-  EXPECT_FLOAT_EQ(1, parameter.states_.back()[1])
-      << "optimal value should be (1, 1)";
+  EXPECT_FLOAT_EQ(0, parameter.states_.front()[3])
+      << "initial value should be (0, 0)";
   EXPECT_FLOAT_EQ(1, parameter.states_.back()[2])
       << "optimal value should be (1, 1)";
+  EXPECT_FLOAT_EQ(1, parameter.states_.back()[3])
+      << "optimal value should be (1, 1)";
   EXPECT_FLOAT_EQ(return_code, 0);
   EXPECT_EQ(19, interrupt.call_count());
 }
diff --git a/src/test/unit/services/optimize/lbfgs_jacobian_test.cpp b/src/test/unit/services/optimize/lbfgs_jacobian_test.cpp
index b6503f10d9a..7ed7365d06f 100644
--- a/src/test/unit/services/optimize/lbfgs_jacobian_test.cpp
+++ b/src/test/unit/services/optimize/lbfgs_jacobian_test.cpp
@@ -35,10 +35,11 @@ TEST_F(ServicesOptimize, with_jacobian) {
   EXPECT_TRUE(logger.find("Optimization terminated normally: "));
   EXPECT_FLOAT_EQ(return_code, 0);
 
-  ASSERT_EQ(2, parameter.names_.size());
+  ASSERT_EQ(3, parameter.names_.size());
   EXPECT_EQ("lp__", parameter.names_[0]);
-  EXPECT_EQ("sigma", parameter.names_[1]);
-  EXPECT_NEAR((3 + std::sqrt(13)) / 2, parameter.states_.back()[1], 0.0001);
+  EXPECT_EQ("converged__", parameter.names_[1]);
+  EXPECT_EQ("sigma", parameter.names_[2]);
+  EXPECT_NEAR((3 + std::sqrt(13)) / 2, parameter.states_.back()[2], 0.0001);
 }
 
 TEST_F(ServicesOptimize, without_jacobian) {
@@ -58,8 +59,9 @@ TEST_F(ServicesOptimize, without_jacobian) {
   EXPECT_TRUE(logger.find("Optimization terminated normally: "));
   EXPECT_FLOAT_EQ(return_code, 0);
 
-  ASSERT_EQ(2, parameter.names_.size());
+  ASSERT_EQ(3, parameter.names_.size());
   EXPECT_EQ("lp__", parameter.names_[0]);
-  EXPECT_EQ("sigma", parameter.names_[1]);
-  EXPECT_NEAR(3, parameter.states_.back()[1], 0.0001);
+  EXPECT_EQ("converged__", parameter.names_[1]);
+  EXPECT_EQ("sigma", parameter.names_[2]);
+  EXPECT_NEAR(3, parameter.states_.back()[2], 0.0001);
 }
diff --git a/src/test/unit/services/optimize/lbfgs_test.cpp b/src/test/unit/services/optimize/lbfgs_test.cpp
index 2a729ec7fc8..243c9b86cc6 100644
--- a/src/test/unit/services/optimize/lbfgs_test.cpp
+++ b/src/test/unit/services/optimize/lbfgs_test.cpp
@@ -40,19 +40,20 @@ TEST_F(ServicesOptimize, rosenbrock) {
 
   EXPECT_EQ("0,0\n", init_ss.str());
 
-  ASSERT_EQ(3, parameter.names_.size());
+  ASSERT_EQ(4, parameter.names_.size());
   EXPECT_EQ("lp__", parameter.names_[0]);
-  EXPECT_EQ("x", parameter.names_[1]);
-  EXPECT_EQ("y", parameter.names_[2]);
+  EXPECT_EQ("converged__", parameter.names_[1]);
+  EXPECT_EQ("x", parameter.names_[2]);
+  EXPECT_EQ("y", parameter.names_[3]);
 
   EXPECT_EQ(23, parameter.states_.size());
-  EXPECT_FLOAT_EQ(0, parameter.states_.front()[1])
-      << "initial value should be (0, 0)";
   EXPECT_FLOAT_EQ(0, parameter.states_.front()[2])
       << "initial value should be (0, 0)";
-  EXPECT_FLOAT_EQ(0.99998301, parameter.states_.back()[1])
+  EXPECT_FLOAT_EQ(0, parameter.states_.front()[3])
+      << "initial value should be (0, 0)";
+  EXPECT_FLOAT_EQ(0.99998301, parameter.states_.back()[2])
       << "optimal value should be (1, 1)";
-  EXPECT_FLOAT_EQ(0.99996597, parameter.states_.back()[2])
+  EXPECT_FLOAT_EQ(0.99996597, parameter.states_.back()[3])
       << "optimal value should be (1, 1)";
   EXPECT_FLOAT_EQ(return_code, 0);
   EXPECT_EQ(22, interrupt.call_count());
diff --git a/src/test/unit/services/optimize/newton_jacobian_test.cpp b/src/test/unit/services/optimize/newton_jacobian_test.cpp
index edd97ff1d3e..00f9de0402c 100644
--- a/src/test/unit/services/optimize/newton_jacobian_test.cpp
+++ b/src/test/unit/services/optimize/newton_jacobian_test.cpp
@@ -32,10 +32,11 @@ TEST_F(ServicesOptimize, withJacobian) {
 
   EXPECT_FLOAT_EQ(return_code, 0);
 
-  ASSERT_EQ(2, parameter.names_.size());
+  ASSERT_EQ(3, parameter.names_.size());
   EXPECT_EQ("lp__", parameter.names_[0]);
-  EXPECT_EQ("sigma", parameter.names_[1]);
-  EXPECT_NEAR((3 + std::sqrt(13)) / 2, parameter.states_.back()[1], 0.001);
+  EXPECT_EQ("converged__", parameter.names_[1]);
+  EXPECT_EQ("sigma", parameter.names_[2]);
+  EXPECT_NEAR((3 + std::sqrt(13)) / 2, parameter.states_.back()[2], 0.001);
   EXPECT_GT(interrupt.call_count(), 0);
 }
 
@@ -54,9 +55,10 @@ TEST_F(ServicesOptimize, withoutJacobian) {
 
   EXPECT_FLOAT_EQ(return_code, 0);
 
-  ASSERT_EQ(2, parameter.names_.size());
+  ASSERT_EQ(3, parameter.names_.size());
   EXPECT_EQ("lp__", parameter.names_[0]);
-  EXPECT_EQ("sigma", parameter.names_[1]);
-  EXPECT_NEAR(3, parameter.states_.back()[1], 0.001);
+  EXPECT_EQ("converged__", parameter.names_[1]);
+  EXPECT_EQ("sigma", parameter.names_[2]);
+  EXPECT_NEAR(3, parameter.states_.back()[2], 0.001);
   EXPECT_GT(interrupt.call_count(), 0);
 }
diff --git a/src/test/unit/services/optimize/newton_test.cpp b/src/test/unit/services/optimize/newton_test.cpp
index daf58006aea..711756922be 100644
--- a/src/test/unit/services/optimize/newton_test.cpp
+++ b/src/test/unit/services/optimize/newton_test.cpp
@@ -36,20 +36,21 @@ TEST_F(ServicesOptimize, rosenbrock) {
   EXPECT_EQ(1, logger.find("Initial log joint probability = -1"));
   EXPECT_EQ(1, logger.find("Iteration  1. Log joint probability ="));
 
-  ASSERT_EQ(3, parameter.names_.size());
+  ASSERT_EQ(4, parameter.names_.size());
   EXPECT_EQ("lp__", parameter.names_[0]);
-  EXPECT_EQ("x", parameter.names_[1]);
-  EXPECT_EQ("y", parameter.names_[2]);
+  EXPECT_EQ("converged__", parameter.names_[1]);
+  EXPECT_EQ("x", parameter.names_[2]);
+  EXPECT_EQ("y", parameter.names_[3]);
 
   EXPECT_GT(parameter.states_.size(), 0);
-  EXPECT_FLOAT_EQ(0, parameter.states_.front()[1])
-      << "initial value should be (0, 0)";
   EXPECT_FLOAT_EQ(0, parameter.states_.front()[2])
       << "initial value should be (0, 0)";
-  EXPECT_NEAR(1, parameter.states_.back()[1], 1e-3)
-      << "optimal value should be (1, 1)";
+  EXPECT_FLOAT_EQ(0, parameter.states_.front()[3])
+      << "initial value should be (0, 0)";
   EXPECT_NEAR(1, parameter.states_.back()[2], 1e-3)
       << "optimal value should be (1, 1)";
+  EXPECT_NEAR(1, parameter.states_.back()[3], 1e-3)
+      << "optimal value should be (1, 1)";
   EXPECT_FLOAT_EQ(return_code, 0);
   EXPECT_LT(0, interrupt.call_count());
 }
@@ -75,16 +76,17 @@ TEST_F(ServicesOptimize, rosenbrock_no_save_iterations) {
 
   EXPECT_EQ("0,0\n", init_ss.str());
 
-  ASSERT_EQ(3, parameter.names_.size());
+  ASSERT_EQ(4, parameter.names_.size());
   EXPECT_EQ("lp__", parameter.names_[0]);
-  EXPECT_EQ("x", parameter.names_[1]);
-  EXPECT_EQ("y", parameter.names_[2]);
+  EXPECT_EQ("converged__", parameter.names_[1]);
+  EXPECT_EQ("x", parameter.names_[2]);
+  EXPECT_EQ("y", parameter.names_[3]);
 
   EXPECT_EQ(1, parameter.states_.size());
-  EXPECT_NEAR(1, parameter.states_.back()[1], 1e-3)
-      << "optimal value should be (1, 1)";
   EXPECT_NEAR(1, parameter.states_.back()[2], 1e-3)
       << "optimal value should be (1, 1)";
+  EXPECT_NEAR(1, parameter.states_.back()[3], 1e-3)
+      << "optimal value should be (1, 1)";
   EXPECT_FLOAT_EQ(return_code, 0);
   EXPECT_LT(0, interrupt.call_count());
 }