diff --git a/Makefile b/Makefile
index 1ba1250..9bf4e44 100644
--- a/Makefile
+++ b/Makefile
@@ -169,7 +169,7 @@ print-%  : ; @echo $* = $($*) ;
 STANC_DL_RETRY = 5
 STANC_DL_DELAY = 10
 STANC3_TEST_BIN_URL ?=
-STANC3_VERSION ?= v2.38.0
+STANC3_VERSION ?= v2.39.0
 
 ifeq ($(OS),Windows_NT)
  OS_TAG := windows
diff --git a/clients/R/R/zzz.R b/clients/R/R/zzz.R
index 90b7400..6da5c62 100644
--- a/clients/R/R/zzz.R
+++ b/clients/R/R/zzz.R
@@ -10,7 +10,7 @@ HMC_SAMPLER_VARIABLES = c(
 
 PATHFINDER_VARIABLES = c("lp_approx__", "lp__", "path__")
 
-OPTIMIZATION_VARIABLES = c("lp__")
+OPTIMIZATION_VARIABLES = c("lp__", "converged__")
 
 LAPLACE_VARIABLES = c("log_p__", "log_g__")
 
diff --git a/clients/julia/src/model.jl b/clients/julia/src/model.jl
index a9a8259..0c431e7 100644
--- a/clients/julia/src/model.jl
+++ b/clients/julia/src/model.jl
@@ -34,7 +34,7 @@ const HMC_SAMPLER_VARIABLES = [
 
 const PATHFINDER_VARIABLES = ["lp_approx__", "lp__", "path__"]
 
-const OPTIMIZE_VARIABLES = ["lp__"]
+const OPTIMIZE_VARIABLES = ["lp__", "converged__"]
 
 const LAPLACE_VARIABLES = ["log_p__", "log_q__"]
 
diff --git a/clients/julia/test/test_laplace.jl b/clients/julia/test/test_laplace.jl
index 88726eb..a9ace56 100644
--- a/clients/julia/test/test_laplace.jl
+++ b/clients/julia/test/test_laplace.jl
@@ -45,7 +45,7 @@
         @testset for jacobian in [true, false]
             out_opt =
                 optimize(simple_jacobian_model; jacobian = jacobian, seed = UInt32(12345))
-            mode_array = out_opt.draws[2:end]
+            mode_array = out_opt.draws[3:end]
             out = laplace_sample(
                 simple_jacobian_model,
                 mode_array;
diff --git a/clients/julia/test/test_optimize.jl b/clients/julia/test/test_optimize.jl
index ad82d55..57c029d 100644
--- a/clients/julia/test/test_optimize.jl
+++ b/clients/julia/test/test_optimize.jl
@@ -96,7 +96,7 @@
 
     @testset "Model without parameters" begin
         out = optimize(empty_model)
-        @test length(out.names) == 1 # lp
+        @test length(out.names) == 2 # lp, converged
     end
 
 
diff --git a/clients/python/tests/test_optimize.py b/clients/python/tests/test_optimize.py
index 4f0a0ec..7501e46 100644
--- a/clients/python/tests/test_optimize.py
+++ b/clients/python/tests/test_optimize.py
@@ -107,7 +107,7 @@ def test_bad_init(bernoulli_model):
 
 def test_model_no_params(empty_model):
     out = empty_model.optimize()
-    assert len(out.parameters) == 1  # lp
+    assert len(out.parameters) == 2  # lp, converged
 
 
 @pytest.mark.parametrize(
diff --git a/clients/python/tinystan/model.py b/clients/python/tinystan/model.py
index a197c88..f41a895 100644
--- a/clients/python/tinystan/model.py
+++ b/clients/python/tinystan/model.py
@@ -66,19 +66,11 @@ def print_callback(msg, size, is_error):
 
 PATHFINDER_VARIABLES = ["lp_approx__", "lp__", "path__"]
 
-OPTIMIZE_VARIABLES = [
-    "lp__",
-]
+OPTIMIZE_VARIABLES = ["lp__", "converged__"]
 
-LAPLACE_VARIABLES = [
-    "log_p__",
-    "log_q__",
-]
+LAPLACE_VARIABLES = ["log_p__", "log_q__"]
 
-FIXED_SAMPLER_VARIABLES = [
-    "lp__",
-    "accept_stat__",
-]
+FIXED_SAMPLER_VARIABLES = ["lp__", "accept_stat__"]
 
 
 class HMCMetric(Enum):
@@ -122,7 +114,7 @@ def preprocess_laplace_inputs(
     if isinstance(mode, StanOutput):
         # handle case of passing optimization output directly
         if len(mode.data.shape) == 1:
-            mode = mode.data[1:]
+            mode = mode.data[len(OPTIMIZE_VARIABLES) :]
         else:
             raise ValueError("Laplace can only be used with Optimization output")
             # mode = mode.create_inits(chains=1, seed=seed)
diff --git a/stan b/stan
index d13c50c..44be14e 160000
--- a/stan
+++ b/stan
@@ -1 +1 @@
-Subproject commit d13c50c0fe7e4a4c387b727b81246cee2541a320
+Subproject commit 44be14e28e8a677719eb821447d37765b3295f7a