diff --git a/src/functions-reference/embedded_laplace.qmd b/src/functions-reference/embedded_laplace.qmd
index 36de2f3e3..5bb8530c0 100644
--- a/src/functions-reference/embedded_laplace.qmd
+++ b/src/functions-reference/embedded_laplace.qmd
@@ -315,6 +315,9 @@ Create a default Laplace options tuple containing `theta_init`.
In `generated quantities`, it is possible to sample from the Laplace
approximation of $p(\theta \mid \phi, y)$ using `laplace_latent_rng`.
+To obtain the approximate mean and Cholesky factor of the covariance
+instead of a random draw, see
+[`laplace_latent_solve`](#return-the-approximate-conditional-mean-and-cholesky-factor).
The signature for `laplace_latent_rng` follows closely
the signature for `laplace_marginal`:
@@ -336,6 +339,55 @@ Samples from the approximate conditional posterior $p(\theta \mid y, \phi)$
and allows the user to tune the control parameters of the approximation.
{{< since 2.39 >}}
+## Return the approximate conditional mean and Cholesky factor
+
+When the parameters of the Gaussian approximation are needed directly,
+rather than a random draw from it, use `laplace_latent_solve`.
+This returns a tuple containing:
+
+1. the posterior mean of the Laplace approximation to $p(\theta \mid y, \phi)$, and
+2. the lower-triangular Cholesky factor of the posterior covariance of that
+ approximation.
+
+The signature for `laplace_latent_solve` follows closely
+the signature for `laplace_marginal`:
+
+
+\index{{\tt \bfseries laplace\_latent\_solve }!{\tt (function likelihood\_function, tuple(...) likelihood\_arguments, int hessian\_block\_size, function covariance\_function, tuple(...) covariance\_arguments): tuple(vector, matrix)}|hyperpage}
+
+`tuple(vector, matrix)` **`laplace_latent_solve`**`(function likelihood_function, tuple(...) likelihood_arguments, int hessian_block_size, function covariance_function, tuple(...) covariance_arguments)`
\newline
+
+Returns the posterior mean and the lower-triangular Cholesky factor of the
+posterior covariance from the Laplace approximation to
+$p(\theta \mid y, \phi)$.
+
+{{< since 2.40 >}}
+
+As with the other embedded Laplace functions, control parameters can be
+specified with a `_tol` variant:
+
+
+\index{{\tt \bfseries laplace\_latent\_solve\_tol }!{\tt (function likelihood\_function, tuple(...) likelihood\_arguments, int hessian\_block\_size, function covariance\_function, tuple(...) covariance\_arguments, tuple(vector, real, int, int, int, int) tolerances): tuple(vector, matrix)}|hyperpage}
+
+`tuple(vector, matrix)` **`laplace_latent_solve_tol`**`(function likelihood_function, tuple(...), int hessian_block_size, function covariance_function, tuple(...), tuple(vector, real, int, int, int, int) tolerances)`
\newline
+
+Returns the posterior mean and the lower-triangular Cholesky factor of the
+posterior covariance from the Laplace approximation to
+$p(\theta \mid y, \phi)$, and allows the user to tune the control parameters
+of the approximation.
+
+{{< since 2.40 >}}
+
+The returned tuple can be unpacked with Stan's tuple indexing, for example:
+
+```stan
+tuple(vector[N], matrix[N, N]) mean_chol
+ = laplace_latent_solve(ll_function, (a, y), hessian_block_size,
+ cov_function, (rho, alpha, x, N, delta));
+vector[N] mu = mean_chol.1;
+matrix[N, N] L = mean_chol.2;
+```
+
## Built-in Laplace marginal likelihood functions
Stan provides convenient wrappers for the embedded Laplace approximation
diff --git a/src/functions-reference/functions_index.qmd b/src/functions-reference/functions_index.qmd
index 4ab04bd2d..952ebe1e7 100644
--- a/src/functions-reference/functions_index.qmd
+++ b/src/functions-reference/functions_index.qmd
@@ -1742,6 +1742,16 @@ pagetitle: Alphabetical Index
- [`(function likelihood_function, tuple(...) likelihood_arguments, int hessian_block_size, function covariance_function, tuple(...) covariance_arguments) : vector`](embedded_laplace.qmd#index-entry-60d4efd35a460871aff97a758f5d366faee26465) (embedded_laplace.html)
+**laplace_latent_solve**:
+
+ - [`(function likelihood_function, tuple(...) likelihood_arguments, int hessian_block_size, function covariance_function, tuple(...) covariance_arguments) : tuple(vector, matrix)`](embedded_laplace.qmd#index-entry-81b0a183fcf12a6986f0736b373ce211aed2fc77) (embedded_laplace.html)
+
+
+**laplace_latent_solve_tol**:
+
+ - [`(function likelihood_function, tuple(...) likelihood_arguments, int hessian_block_size, function covariance_function, tuple(...) covariance_arguments, tuple(vector, real, int, int, int, int) tolerances) : tuple(vector, matrix)`](embedded_laplace.qmd#index-entry-1dcf9701db2ac630332a16b39610117eff866f6a) (embedded_laplace.html)
+
+
**laplace_latent_tol_bernoulli_logit_rng**:
- [`(array[] int y, array[] int y_index, vector m, data int hessian_block_size, function covariance_function, tuple(...) covariance_arguments, tuple(vector, real, int, int, int, int) tolerances) : vector`](embedded_laplace.qmd#index-entry-4566ad7c6cd5f0ee057539a6e2a0bea28554f082) (embedded_laplace.html)