diff --git a/DESCRIPTION b/DESCRIPTION
index e6f376a4..1c60c03c 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -16,7 +16,8 @@ Authors@R:
       person(given = c("William", "Michael"), family = "Landau", role = "ctb",
            email = "will.landau@gmail.com", comment = c(ORCID = "0000-0003-1878-3253")),
       person(given = "Jacob", family = "Socolar", role = "ctb"), 
-      person(given = "Martin", family = "Modrák", role = "ctb"))
+      person(given = "Martin", family = "Modrák", role = "ctb"),
+      person(given = "Steve", family = "Bronder", role = "ctb"))
 Description: A lightweight interface to 'Stan' <https://mc-stan.org>.
     The 'CmdStanR' interface is an alternative to 'RStan' that calls the command
     line interface for compilation and running algorithms instead of interfacing
diff --git a/NAMESPACE b/NAMESPACE
index 4ab6a168..c8a6217d 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -4,6 +4,7 @@ S3method(as_draws,CmdStanGQ)
 S3method(as_draws,CmdStanLaplace)
 S3method(as_draws,CmdStanMCMC)
 S3method(as_draws,CmdStanMLE)
+S3method(as_draws,CmdStanPathfinder)
 S3method(as_draws,CmdStanVB)
 export(as_cmdstan_fit)
 export(as_draws)
diff --git a/R/args.R b/R/args.R
index 4ac5ab76..123022a5 100644
--- a/R/args.R
+++ b/R/args.R
@@ -16,6 +16,7 @@
 #' * `OptimizeArgs`: stores arguments specific to `method=optimize`.
 #' * `LaplaceArgs`: stores arguments specific to `method=laplace`.
 #' * `VariationalArgs`: stores arguments specific to `method=variational`
+#' * `PathfinderArgs`: stores arguments specific to `method=pathfinder`
 #' * `GenerateQuantitiesArgs`: stores arguments specific to `method=generate_quantities`
 #' * `DiagnoseArgs`: stores arguments specific to `method=diagnose`
 #'
@@ -41,7 +42,8 @@ CmdStanArgs <- R6::R6Class(
                           output_basename = NULL,
                           sig_figs = NULL,
                           opencl_ids = NULL,
-                          model_variables = NULL) {
+                          model_variables = NULL,
+                          num_threads = NULL) {
 
       self$model_name <- model_name
       self$stan_code <- stan_code
@@ -82,6 +84,7 @@ CmdStanArgs <- R6::R6Class(
       }
       self$init <- init
       self$opencl_ids <- opencl_ids
+      self$num_threads = NULL
       self$method_args$validate(num_procs = length(self$proc_ids))
       self$validate()
     },
@@ -179,6 +182,9 @@ CmdStanArgs <- R6::R6Class(
       if (!is.null(self$opencl_ids)) {
         args$opencl <- c("opencl", paste0("platform=", self$opencl_ids[1]), paste0("device=", self$opencl_ids[2]))
       }
+      if (!is.null(self$num_threads)) {
+        num_threads <- c(args$output, paste0("num_threads=", self$num_threads))
+      }
       args <- do.call(c, append(args, list(use.names = FALSE)))
       self$method_args$compose(idx, args)
     },
@@ -541,6 +547,75 @@ VariationalArgs <- R6::R6Class(
   )
 )
 
+# PathfinderArgs ---------------------------------------------------------
+
+PathfinderArgs <- R6::R6Class(
+  "PathfinderArgs",
+  lock_objects = FALSE,
+  public = list(
+    method = "pathfinder",
+      initialize = function(init_alpha = NULL,
+                            tol_obj = NULL,
+                            tol_rel_obj = NULL,
+                            tol_grad = NULL,
+                            tol_rel_grad = NULL,
+                            tol_param = NULL,
+                            history_size = NULL,
+                            single_path_draws = NULL,
+                            draws = NULL,
+                            num_paths = NULL,
+                            max_lbfgs_iters = NULL,
+                            num_elbo_draws = NULL,
+                            save_single_paths = NULL) {
+        self$init_alpha <- init_alpha
+        self$tol_obj <- tol_obj
+        self$tol_rel_obj <- tol_rel_obj
+        self$tol_grad <- tol_grad
+        self$tol_rel_grad <- tol_rel_grad
+        self$tol_param <- tol_param
+        self$history_size <- history_size
+        self$num_psis_draws <- draws
+        self$num_draws <- single_path_draws
+        self$num_paths <- num_paths
+        self$max_lbfgs_iters <- max_lbfgs_iters
+        self$num_elbo_draws <- num_elbo_draws
+        self$save_single_paths <- save_single_paths
+      invisible(self)
+    },
+
+    validate = function(num_procs) {
+      validate_pathfinder_args(self)
+    },
+
+    # Compose arguments to CmdStan command for pathfinder-specific
+    # non-default arguments. Works the same way as compose for sampler args,
+    # but `idx` (multiple pathfinders are handled in cmdstan)
+    compose = function(idx = NULL, args = NULL) {
+      .make_arg <- function(arg_name) {
+        compose_arg(self, arg_name, idx = NULL)
+      }
+        new_args <- list(
+          "method=pathfinder",
+          .make_arg("init_alpha"),
+          .make_arg("tol_obj"),
+          .make_arg("tol_rel_obj"),
+          .make_arg("tol_grad"),
+          .make_arg("tol_rel_grad"),
+          .make_arg("tol_param"),
+          .make_arg("history_size"),
+          .make_arg("num_psis_draws"),
+          .make_arg("num_draws"),
+          .make_arg("num_paths"),
+          .make_arg("max_lbfgs_iters"),
+          .make_arg("num_elbo_draws"),
+          .make_arg("save_single_paths")
+        )
+      new_args <- do.call(c, new_args)
+      c(args, new_args)
+    }
+  )
+)
+
 # DiagnoseArgs -------------------------------------------------------------
 
 DiagnoseArgs <- R6::R6Class(
@@ -854,6 +929,59 @@ validate_variational_args <- function(self) {
   invisible(TRUE)
 }
 
+#' Validate arguments for pathfinder inference
+#' @noRd
+#' @param self A `PathfinderArgs` object.
+#' @return `TRUE` invisibly unless an error is thrown.
+validate_pathfinder_args <- function(self) {
+
+  checkmate::assert_integerish(self$max_lbfgs_iters, lower = 1, null.ok = TRUE, len = 1)
+  if (!is.null(self$max_lbfgs_iters)) {
+    self$iter <- as.integer(self$max_lbfgs_iters)
+  }
+  checkmate::assert_integerish(self$num_paths, lower = 1, null.ok = TRUE,
+                               len = 1)
+  if (!is.null(self$num_paths)) {
+    self$num_paths <- as.integer(self$num_paths)
+  }
+  checkmate::assert_integerish(self$num_draws, lower = 1, null.ok = TRUE,
+                               len = 1, .var.name = "single_path_draws")
+  if (!is.null(self$num_draws)) {
+    self$num_draws <- as.integer(self$num_draws)
+  }
+  checkmate::assert_integerish(self$num_psis_draws, lower = 1, null.ok = TRUE,
+                               len = 1, .var.name = "draws")
+  if (!is.null(self$num_psis_draws)) {
+    self$num_psis_draws <- as.integer(self$num_psis_draws)
+  }
+  checkmate::assert_integerish(self$num_elbo_draws, lower = 1, null.ok = TRUE, len = 1)
+  if (!is.null(self$num_elbo_draws)) {
+    self$num_elbo_draws <- as.integer(self$num_elbo_draws)
+  }
+  if (!is.null(self$save_single_paths) && is.logical(self$save_single_paths)) {
+    self$save_single_paths = as.integer(self$save_single_paths)
+  }
+  checkmate::assert_integerish(self$save_single_paths, null.ok = TRUE,
+                               lower = 0, upper = 1, len = 1)
+  if (!is.null(self$save_single_paths)) {
+    self$save_single_paths <- 0
+  }
+
+
+  # check args only available for lbfgs and bfgs
+  bfgs_args <- c("init_alpha", "tol_obj", "tol_rel_obj", "tol_grad", "tol_rel_grad", "tol_param")
+  for (arg in bfgs_args) {
+    checkmate::assert_number(self[[arg]], .var.name = arg, lower = 0, null.ok = TRUE)
+  }
+
+  if (!is.null(self$history_size)) {
+      checkmate::assert_integerish(self$history_size, lower = 1, len = 1, null.ok = FALSE)
+      self$history_size <- as.integer(self$history_size)
+  }
+
+  invisible(TRUE)
+}
+
 
 # Validation helpers ------------------------------------------------------
 
diff --git a/R/csv.R b/R/csv.R
index 33095850..f2e674f8 100644
--- a/R/csv.R
+++ b/R/csv.R
@@ -280,6 +280,10 @@ read_cmdstan_csv <- function(files,
     if (length(variables) > 0) {
       draws_list_id <- length(draws) + 1
       warmup_draws_list_id <- length(warmup_draws) + 1
+      if (metadata$method == "pathfinder") {
+        metadata$variables = union(metadata$sampler_diagnostics, metadata$variables)
+        variables = union(metadata$sampler_diagnostics, variables)
+      }
       suppressWarnings(
         draws[[draws_list_id]] <- data.table::fread(
           cmd = fread_cmd,
@@ -445,6 +449,21 @@ read_cmdstan_csv <- function(files,
       metadata = metadata,
       generated_quantities = draws
     )
+  } else if (metadata$method == "pathfinder") {
+    if (is.null(format)) {
+      format <- "draws_matrix"
+    }
+    as_draws_format <- as_draws_format_fun(format)
+    if (length(draws) == 0) {
+      pathfinder_draws <- NULL
+    } else {
+      pathfinder_draws <- do.call(as_draws_format, list(draws[[1]][, colnames(draws[[1]]), drop = FALSE]))
+      posterior::variables(pathfinder_draws) <- repaired_variables
+    }
+    list(
+      metadata = metadata,
+      draws = pathfinder_draws
+    )
   }
 }
 
@@ -477,6 +496,7 @@ as_cmdstan_fit <- function(files, check_diagnostics = TRUE, format = getOption("
     "sample" = CmdStanMCMC_CSV$new(csv_contents, files, check_diagnostics),
     "optimize" = CmdStanMLE_CSV$new(csv_contents, files),
     "variational" = CmdStanVB_CSV$new(csv_contents, files),
+    "pathfinder" = CmdStanPathfinder_CSV$new(csv_contents, files),
     "laplace" = CmdStanLaplace_CSV$new(csv_contents, files)
   )
 }
@@ -576,6 +596,23 @@ CmdStanVB_CSV <- R6::R6Class(
   private = list(output_files_ = NULL)
 )
 
+CmdStanPathfinder_CSV <- R6::R6Class(
+  classname = "CmdStanPathfinder_CSV",
+  inherit = CmdStanPathfinder,
+  public = list(
+    initialize = function(csv_contents, files) {
+      private$output_files_ <- files
+      private$draws_ <- csv_contents$draws
+      private$metadata_ <- csv_contents$metadata
+    },
+    output_files = function(...) {
+      private$output_files_
+    }
+  ),
+  private = list(output_files_ = NULL)
+)
+
+
 # these methods are unavailable because there's no CmdStanRun object
 unavailable_methods_CmdStanFit_CSV <- c(
     "cmdstan_diagnose", "cmdstan_summary",
@@ -642,7 +679,7 @@ read_csv_metadata <- function(csv_file) {
       "\""
     )
   } else {
-    fread_cmd <- paste0("grep '^[#a-zA-Z]' --color=never '", csv_file, "'")
+    fread_cmd <- paste0("grep '^[#a-zA-Z]' --color=never '", path.expand(csv_file), "'")
   }
   suppressWarnings(
     metadata <- data.table::fread(
diff --git a/R/example.R b/R/example.R
index 924c961f..5bb98ae3 100644
--- a/R/example.R
+++ b/R/example.R
@@ -50,7 +50,7 @@
 #'
 cmdstanr_example <-
   function(example = c("logistic", "schools", "schools_ncp"),
-           method = c("sample", "optimize", "laplace", "variational", "diagnose"),
+           method = c("sample", "optimize", "laplace", "variational", "pathfinder", "diagnose"),
            ...,
            quiet = TRUE,
            force_recompile = getOption("cmdstanr_force_recompile", default = FALSE)) {
diff --git a/R/fit.R b/R/fit.R
index 71aa889e..164e420e 100644
--- a/R/fit.R
+++ b/R/fit.R
@@ -2024,6 +2024,80 @@ CmdStanVB <- R6::R6Class(
 )
 CmdStanVB$set("public", name = "lp_approx", value = lp_approx)
 
+# CmdStanPathfinder ---------------------------------------------------------------
+#' CmdStanPathfinder objects
+#'
+#' @name CmdStanPathfinder
+#' @family fitted model objects
+#' @template seealso-docs
+#'
+#' @description A `CmdStanPathfinder` object is the fitted model object returned by the
+#'   [`$pathfinder()`][model-method-pathfinder] method of a
+#'   [`CmdStanModel`] object.
+#'
+#' @section Methods: `CmdStanPathfinder` objects have the following associated methods,
+#'   all of which have their own (linked) documentation pages.
+#'
+#'  ## Extract contents of fitted model object
+#'
+#'  |**Method**|**Description**|
+#'  |:----------|:---------------|
+#'  [`$draws()`][fit-method-draws]  |  Return approximate posterior draws as a [`draws_matrix`][posterior::draws_matrix]. |
+#'  [`$lp()`][fit-method-lp]  |  Return the total log probability density (`target`) computed in the model block of the Stan program. |
+#'  [`$lp_approx()`][fit-method-lp]  |  Return the log density of the approximation to the posterior. |
+#'  [`$init()`][fit-method-init] |  Return user-specified initial values. |
+#'  [`$metadata()`][fit-method-metadata] | Return a list of metadata gathered from the CmdStan CSV files. |
+#'  [`$code()`][fit-method-code] | Return Stan code as a character vector. |
+#'
+#'  ## Summarize inferences
+#'
+#'  |**Method**|**Description**|
+#'  |:----------|:---------------|
+#'  [`$summary()`][fit-method-summary]  | Run [`posterior::summarise_draws()`][posterior::draws_summary]. |
+#'  [`$cmdstan_summary()`][fit-method-cmdstan_summary] |  Run and print CmdStan's `bin/stansummary`. |
+#'
+#'  ## Save fitted model object and temporary files
+#'
+#'  |**Method**|**Description**|
+#'  |:----------|:---------------|
+#'  [`$save_object()`][fit-method-save_object] |  Save fitted model object to a file. |
+#'  [`$save_output_files()`][fit-method-save_output_files] |  Save output CSV files to a specified location. |
+#'  [`$save_data_file()`][fit-method-save_data_file] |  Save JSON data file to a specified location. |
+#'  [`$save_latent_dynamics_files()`][fit-method-save_latent_dynamics_files] |  Save diagnostic CSV files to a specified location. |
+#'
+#'  ## Report run times, console output, return codes
+#'
+#'  |**Method**|**Description**|
+#'  |:----------|:---------------|
+#'  [`$time()`][fit-method-time]  |  Report the total run time. |
+#'  [`$output()`][fit-method-output]  |  Pretty print the output that was printed to the console. |
+#'  [`$return_codes()`][fit-method-return_codes]  |  Return the return codes from the CmdStan runs. |
+#'
+CmdStanPathfinder <- R6::R6Class(
+  classname = "CmdStanPathfinder",
+  inherit = CmdStanFit,
+  public = list(),
+  private = list(
+    # inherits draws_ and metadata_ slots from CmdStanFit
+    read_csv_ = function(format = getOption("cmdstanr_draws_format", "draws_matrix")) {
+      if (!length(self$output_files(include_failed = FALSE))) {
+        stop("Pathfinder failed. Unable to retrieve the draws.", call. = FALSE)
+      }
+      csv_contents <- read_cmdstan_csv(self$output_files(), format = format)
+      private$draws_ <- csv_contents$draws
+      private$metadata_ <- csv_contents$metadata
+      invisible(self)
+    }
+  )
+)
+
+#' @rdname fit-method-lp
+lp_approx <- function() {
+  as.numeric(self$draws()[, "lp_approx__"])
+}
+CmdStanPathfinder$set("public", name = "lp_approx", value = lp_approx)
+
+
 
 # CmdStanGQ ---------------------------------------------------------------
 #' CmdStanGQ objects
@@ -2290,3 +2364,9 @@ as_draws.CmdStanVB <- function(x, ...) {
 as_draws.CmdStanGQ <- function(x, ...) {
   x$draws(...)
 }
+
+#' @rdname as_draws.CmdStanMCMC
+#' @export
+as_draws.CmdStanPathfinder <- function(x, ...) {
+  x$draws(...)
+}
diff --git a/R/model.R b/R/model.R
index c8d6ea60..f9a0e74d 100644
--- a/R/model.R
+++ b/R/model.R
@@ -99,13 +99,15 @@
 #' fit_laplace <- mod$laplace(data = my_data_file, mode = fit_optim, draws = 2000)
 #' fit_laplace$summary()
 #'
-#' # Run 'variational' method to approximate the posterior (default is meanfield ADVI)
+#' # Run 'variational' method to use ADVI to approximate posterior
 #' fit_vb <- mod$variational(data = stan_data, seed = 123)
 #' fit_vb$summary()
-#'
-#' # Plot approximate posterior using bayesplot
 #' mcmc_hist(fit_vb$draws("theta"))
 #'
+#' # Run 'pathfinder' method, a new alternative to the variational method
+#' fit_pf <- mod$pathfinder(data = stan_data, seed = 123)
+#' fit_pf$summary()
+#' mcmc_hist(fit_pf$draws("theta"))
 #'
 #' # Specifying initial values as a function
 #' fit_mcmc_w_init_fun <- mod$sample(
@@ -205,6 +207,7 @@ cmdstan_model <- function(stan_file = NULL, exe_file = NULL, compile = TRUE, ...
 #'  [`$sample_mpi()`][model-method-sample_mpi] |  Run CmdStan's `"sample"` method with [MPI](https://mc-stan.org/math/mpi.html), return [`CmdStanMCMC`] object. |
 #'  [`$optimize()`][model-method-optimize] |  Run CmdStan's `"optimize"` method, return [`CmdStanMLE`] object. |
 #'  [`$variational()`][model-method-variational] |  Run CmdStan's `"variational"` method, return [`CmdStanVB`] object. |
+#'  [`$pathfinder()`][model-method-pathfinder] |  Run CmdStan's `"pathfinder"` method, return [`CmdStanPathfinder`] object. |
 #'  [`$generate_quantities()`][model-method-generate-quantities] |  Run CmdStan's `"generate quantities"` method, return [`CmdStanGQ`] object. |
 #'
 #' @template seealso-docs
@@ -1426,8 +1429,8 @@ CmdStanModel$set("public", name = "sample_mpi", value = sample_mpi)
 #' @description The `$optimize()` method of a [`CmdStanModel`] object runs
 #'   Stan's optimizer to obtain a (penalized) maximum likelihood estimate or a
 #'   maximum a posteriori estimate (if `jacobian=TRUE`). See the
-#'   [Maximum Likelihood Estimation](https://mc-stan.org/docs/cmdstan-guide/maximum-likelihood-estimation.html)
-#'   section of the CmdStan User's Guide for more details.
+#'   [CmdStan User's Guide](https://mc-stan.org/docs/cmdstan-guide/index.html)
+#'   for more details.
 #'
 #'   Any argument left as `NULL` will default to the default value used by the
 #'   installed version of CmdStan. See the [CmdStan User’s
@@ -1550,13 +1553,11 @@ CmdStanModel$set("public", name = "optimize", value = optimize)
 #'   likelihood estimate (MLE), the sample provides an estimate of the standard
 #'   error of the likelihood. Whether the mode is the MAP or MLE depends on
 #'   the value of the `jacobian` argument when running optimization. See the
-#'   [Laplace Sampling](https://mc-stan.org/docs/cmdstan-guide/laplace-sampling.html)
-#'   section of the CmdStan User's Guide for more details.
+#'   [CmdStan User’s Guide](https://mc-stan.org/docs/cmdstan-guide/)
+#'   for more details.
 #'
 #'   Any argument left as `NULL` will default to the default value used by the
-#'   installed version of CmdStan. See the
-#'   [CmdStan User’s Guide](https://mc-stan.org/docs/cmdstan-guide/)
-#'   for more details on the default arguments.
+#'   installed version of CmdStan.
 #'
 #' @template model-common-args
 #' @inheritParams model-method-optimize
@@ -1707,21 +1708,17 @@ CmdStanModel$set("public", name = "laplace", value = laplace)
 #' @family CmdStanModel methods
 #'
 #' @description The `$variational()` method of a [`CmdStanModel`] object runs
-#'   Stan's variational Bayes (ADVI) algorithms.
-#'
-#'   Any argument left as `NULL` will default to the default value used by the
-#'   installed version of CmdStan. See the
+#'   Stan's Automatic Differentiation Variational Inference (ADVI) algorithms.
+#'   The approximation is a Gaussian in the unconstrained variable space. Stan
+#'   implements two ADVI algorithms: the `algorithm="meanfield"` option uses a
+#'   fully factorized Gaussian for the approximation; the `algorithm="fullrank"`
+#'   option uses a Gaussian with a full-rank covariance matrix for the
+#'   approximation. See the
 #'   [CmdStan User’s Guide](https://mc-stan.org/docs/cmdstan-guide/)
 #'   for more details.
 #'
-#' @details CmdStan can fit a variational approximation to the posterior. The
-#'   approximation is a Gaussian in the unconstrained variable space. Stan
-#'   implements two variational algorithms. The `algorithm="meanfield"` option
-#'   uses a fully factorized Gaussian for the approximation. The
-#'   `algorithm="fullrank"` option uses a Gaussian with a full-rank covariance
-#'   matrix for the approximation.
-#'
-#'   -- [*CmdStan Interface User's Guide*](https://github.com/stan-dev/cmdstan/releases/latest)
+#'   Any argument left as `NULL` will default to the default value used by the
+#'   installed version of CmdStan.
 #'
 #' @template model-common-args
 #' @param threads (positive integer) If the model was
@@ -1821,6 +1818,131 @@ variational <- function(data = NULL,
 }
 CmdStanModel$set("public", name = "variational", value = variational)
 
+
+#' Run Stan's Pathfinder Variational Inference Algorithm
+#'
+#' @name model-method-pathfinder
+#' @aliases pathfinder
+#' @family CmdStanModel methods
+#'
+#' @description The `$pathfinder()` method of a [`CmdStanModel`] object runs
+#'   Stan's Pathfinder algorithms. Pathfinder is a variational method for
+#'   approximately sampling from differentiable log densities. Starting from a
+#'   random initialization, Pathfinder locates normal approximations to the
+#'   target density along a quasi-Newton optimization path, with local
+#'   covariance estimated using the negative inverse Hessian estimates produced
+#'   by the L-BFGS optimizer. Pathfinder returns draws from the Gaussian
+#'   approximation with the lowest estimated Kullback-Leibler (KL) divergence to
+#'   the true posterior. See the
+#'   [CmdStan User’s Guide](https://mc-stan.org/docs/cmdstan-guide/)
+#'   for more details.
+#'
+#'   Any argument left as `NULL` will default to the default value used by the
+#'   installed version of CmdStan
+#'
+#' @template model-common-args
+#' @param num_threads (positive integer) If the model was
+#'   [compiled][model-method-compile] with threading support, the number of
+#'   threads to use in parallelized sections (e.g., for multi-path pathfinder
+#'   as well as `reduce_sum`).
+#' @param init_alpha (positive real) The initial step size parameter.
+#' @param tol_obj (positive real) Convergence tolerance on changes in objective function value.
+#' @param tol_rel_obj (positive real) Convergence tolerance on relative changes in objective function value.
+#' @param tol_grad (positive real) Convergence tolerance on the norm of the gradient.
+#' @param tol_rel_grad (positive real) Convergence tolerance on the relative norm of the gradient.
+#' @param tol_param (positive real) Convergence tolerance on changes in parameter value.
+#' @param history_size (positive integer) The size of the history used when
+#'   approximating the Hessian.
+#' @param single_path_draws (positive integer) Number of draws a single
+#'   pathfinder should return. The number of draws PSIS sampling samples from
+#'   will be equal to `single_path_draws * num_paths`.
+#' @param draws (positive integer) Number of draws to return after performing
+#'   pareto smooted importance sampling (PSIS). This must be smaller than
+#'   `single_path_draws * num_paths`.
+#' @param num_paths (positive integer) Number of single pathfinders to run.
+#' @param max_lbfgs_iters (positive integer) The maximum number of iterations
+#'   for LBFGS.
+#' @param num_elbo_draws (positive integer) Number of draws to make when
+#'   calculating the ELBO of the approximation at each iteration of LBFGS.
+#' @param save_single_paths (logical) Whether to save the results of single
+#'   pathfinder runs in multi-pathfinder.
+#' @return A [`CmdStanPathfinder`] object.
+#'
+#' @template seealso-docs
+#' @inherit cmdstan_model examples
+#'
+pathfinder <- function(data = NULL,
+                       seed = NULL,
+                       refresh = NULL,
+                       init = NULL,
+                       save_latent_dynamics = FALSE,
+                       output_dir = NULL,
+                       output_basename = NULL,
+                       sig_figs = NULL,
+                       opencl_ids = NULL,
+                       num_threads = NULL,
+                       init_alpha = NULL,
+                       tol_obj = NULL,
+                       tol_rel_obj = NULL,
+                       tol_grad = NULL,
+                       tol_rel_grad = NULL,
+                       tol_param = NULL,
+                       history_size = NULL,
+                       single_path_draws = NULL,
+                       draws = NULL,
+                       num_paths = NULL,
+                       max_lbfgs_iters = NULL,
+                       num_elbo_draws = NULL,
+                       save_single_paths = NULL) {
+  procs <- CmdStanProcs$new(
+    num_procs = 1,
+    show_stdout_messages = (is.null(refresh) || refresh != 0),
+    threads_per_proc = assert_valid_threads(num_threads, self$cpp_options())
+  )
+  model_variables <- NULL
+  if (is_variables_method_supported(self)) {
+    model_variables <- self$variables()
+  }
+  pathfinder_args <- PathfinderArgs$new(
+   init_alpha = init_alpha,
+   tol_obj = tol_obj,
+   tol_rel_obj = tol_rel_obj,
+   tol_grad = tol_grad,
+   tol_rel_grad = tol_rel_grad,
+   tol_param = tol_param,
+   history_size = history_size,
+   draws = draws,
+   single_path_draws = single_path_draws,
+   num_paths = num_paths,
+   max_lbfgs_iters = max_lbfgs_iters,
+   num_elbo_draws = num_elbo_draws,
+   save_single_paths = save_single_paths)
+  args <- CmdStanArgs$new(
+    method_args = pathfinder_args,
+    stan_file = self$stan_file(),
+    stan_code = suppressWarnings(self$code()),
+    model_name = self$model_name(),
+    exe_file = self$exe_file(),
+    proc_ids = 1,
+    data_file = process_data(data, model_variables),
+    save_latent_dynamics = save_latent_dynamics,
+    seed = seed,
+    init = init,
+    refresh = refresh,
+    output_dir = output_dir,
+    output_basename = output_basename,
+    sig_figs = sig_figs,
+    opencl_ids = assert_valid_opencl(opencl_ids, self$cpp_options()),
+    model_variables = model_variables,
+    num_threads = num_threads
+  )
+  runset <- CmdStanRun$new(args, procs)
+  runset$run_cmdstan()
+  CmdStanPathfinder$new(runset)
+}
+CmdStanModel$set("public", name = "pathfinder", value = pathfinder)
+
+
 #' Run Stan's standalone generated quantities method
 #'
 #' @name model-method-generate-quantities
diff --git a/R/run.R b/R/run.R
index b674a05e..16a0c97b 100644
--- a/R/run.R
+++ b/R/run.R
@@ -262,7 +262,7 @@ CmdStanRun <- R6::R6Class(
     },
 
     time = function() {
-      if (self$method() %in% c("laplace", "optimize", "variational")) {
+      if (self$method() %in% c("laplace", "optimize", "variational", "pathfinder")) {
         time <- list(total = self$procs$total_time())
       } else if (self$method() == "generate_quantities") {
         chain_time <- data.frame(
@@ -518,6 +518,10 @@ CmdStanRun$set("private", name = "run_generate_quantities_", value = .run_genera
     if (procs$proc_state(id = id) > 3) {
       successful_fit <- TRUE
     }
+  } else if (self$method() == "pathfinder") {
+    if (procs$proc_state(id = id) > 3 | procs$get_proc(id)$get_exit_status() == 0) {
+      successful_fit <- TRUE
+    }
   } else if (procs$get_proc(id)$get_exit_status() == 0) {
     successful_fit <- TRUE
   }
@@ -533,6 +537,7 @@ CmdStanRun$set("private", name = "run_generate_quantities_", value = .run_genera
 CmdStanRun$set("private", name = "run_optimize_", value = .run_other)
 CmdStanRun$set("private", name = "run_laplace_", value = .run_other)
 CmdStanRun$set("private", name = "run_variational_", value = .run_other)
+CmdStanRun$set("private", name = "run_pathfinder_", value = .run_other)
 
 .run_diagnose <- function() {
   procs <- self$procs
diff --git a/man/CmdStanDiagnose.Rd b/man/CmdStanDiagnose.Rd
index 70930db1..8c786e88 100644
--- a/man/CmdStanDiagnose.Rd
+++ b/man/CmdStanDiagnose.Rd
@@ -45,6 +45,7 @@ Other fitted model objects:
 \code{\link{CmdStanLaplace}},
 \code{\link{CmdStanMCMC}},
 \code{\link{CmdStanMLE}},
+\code{\link{CmdStanPathfinder}},
 \code{\link{CmdStanVB}}
 }
 \concept{fitted model objects}
diff --git a/man/CmdStanGQ.Rd b/man/CmdStanGQ.Rd
index 6ebe2d9c..73fab133 100644
--- a/man/CmdStanGQ.Rd
+++ b/man/CmdStanGQ.Rd
@@ -106,6 +106,7 @@ Other fitted model objects:
 \code{\link{CmdStanLaplace}},
 \code{\link{CmdStanMCMC}},
 \code{\link{CmdStanMLE}},
+\code{\link{CmdStanPathfinder}},
 \code{\link{CmdStanVB}}
 }
 \concept{fitted model objects}
diff --git a/man/CmdStanLaplace.Rd b/man/CmdStanLaplace.Rd
index 50d782a6..4d16723e 100644
--- a/man/CmdStanLaplace.Rd
+++ b/man/CmdStanLaplace.Rd
@@ -67,6 +67,7 @@ Other fitted model objects:
 \code{\link{CmdStanGQ}},
 \code{\link{CmdStanMCMC}},
 \code{\link{CmdStanMLE}},
+\code{\link{CmdStanPathfinder}},
 \code{\link{CmdStanVB}}
 }
 \concept{fitted model objects}
diff --git a/man/CmdStanMCMC.Rd b/man/CmdStanMCMC.Rd
index c74ff41f..2148b7a6 100644
--- a/man/CmdStanMCMC.Rd
+++ b/man/CmdStanMCMC.Rd
@@ -89,6 +89,7 @@ Other fitted model objects:
 \code{\link{CmdStanGQ}},
 \code{\link{CmdStanLaplace}},
 \code{\link{CmdStanMLE}},
+\code{\link{CmdStanPathfinder}},
 \code{\link{CmdStanVB}}
 }
 \concept{fitted model objects}
diff --git a/man/CmdStanMLE.Rd b/man/CmdStanMLE.Rd
index 141e45aa..82c11910 100644
--- a/man/CmdStanMLE.Rd
+++ b/man/CmdStanMLE.Rd
@@ -79,6 +79,7 @@ Other fitted model objects:
 \code{\link{CmdStanGQ}},
 \code{\link{CmdStanLaplace}},
 \code{\link{CmdStanMCMC}},
+\code{\link{CmdStanPathfinder}},
 \code{\link{CmdStanVB}}
 }
 \concept{fitted model objects}
diff --git a/man/CmdStanModel.Rd b/man/CmdStanModel.Rd
index 2644784e..7f2c1f9f 100644
--- a/man/CmdStanModel.Rd
+++ b/man/CmdStanModel.Rd
@@ -47,6 +47,7 @@ methods, many of which have their own (linked) documentation pages:
    \code{\link[=model-method-sample_mpi]{$sample_mpi()}} \tab Run CmdStan's \code{"sample"} method with \href{https://mc-stan.org/math/mpi.html}{MPI}, return \code{\link{CmdStanMCMC}} object. \cr
    \code{\link[=model-method-optimize]{$optimize()}} \tab Run CmdStan's \code{"optimize"} method, return \code{\link{CmdStanMLE}} object. \cr
    \code{\link[=model-method-variational]{$variational()}} \tab Run CmdStan's \code{"variational"} method, return \code{\link{CmdStanVB}} object. \cr
+   \code{\link[=model-method-pathfinder]{$pathfinder()}} \tab Run CmdStan's \code{"pathfinder"} method, return \code{\link{CmdStanPathfinder}} object. \cr
    \code{\link[=model-method-generate-quantities]{$generate_quantities()}} \tab Run CmdStan's \code{"generate quantities"} method, return \code{\link{CmdStanGQ}} object. \cr
 }
 
@@ -117,13 +118,15 @@ fit_optim <- mod$optimize(data = my_data_file, jacobian = TRUE)
 fit_laplace <- mod$laplace(data = my_data_file, mode = fit_optim, draws = 2000)
 fit_laplace$summary()
 
-# Run 'variational' method to approximate the posterior (default is meanfield ADVI)
+# Run 'variational' method to use ADVI to approximate posterior
 fit_vb <- mod$variational(data = stan_data, seed = 123)
 fit_vb$summary()
-
-# Plot approximate posterior using bayesplot
 mcmc_hist(fit_vb$draws("theta"))
 
+# Run 'pathfinder' method, a new alternative to the variational method
+fit_pf <- mod$pathfinder(data = stan_data, seed = 123)
+fit_pf$summary()
+mcmc_hist(fit_pf$draws("theta"))
 
 # Specifying initial values as a function
 fit_mcmc_w_init_fun <- mod$sample(
diff --git a/man/CmdStanPathfinder.Rd b/man/CmdStanPathfinder.Rd
new file mode 100644
index 00000000..970174dc
--- /dev/null
+++ b/man/CmdStanPathfinder.Rd
@@ -0,0 +1,73 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/fit.R
+\name{CmdStanPathfinder}
+\alias{CmdStanPathfinder}
+\title{CmdStanPathfinder objects}
+\description{
+A \code{CmdStanPathfinder} object is the fitted model object returned by the
+\code{\link[=model-method-pathfinder]{$pathfinder()}} method of a
+\code{\link{CmdStanModel}} object.
+}
+\section{Methods}{
+ \code{CmdStanPathfinder} objects have the following associated methods,
+all of which have their own (linked) documentation pages.
+\subsection{Extract contents of fitted model object}{\tabular{ll}{
+   \strong{Method} \tab \strong{Description} \cr
+   \code{\link[=fit-method-draws]{$draws()}} \tab Return approximate posterior draws as a \code{\link[posterior:draws_matrix]{draws_matrix}}. \cr
+   \code{\link[=fit-method-lp]{$lp()}} \tab Return the total log probability density (\code{target}) computed in the model block of the Stan program. \cr
+   \code{\link[=fit-method-lp]{$lp_approx()}} \tab Return the log density of the approximation to the posterior. \cr
+   \code{\link[=fit-method-init]{$init()}} \tab Return user-specified initial values. \cr
+   \code{\link[=fit-method-metadata]{$metadata()}} \tab Return a list of metadata gathered from the CmdStan CSV files. \cr
+   \code{\link[=fit-method-code]{$code()}} \tab Return Stan code as a character vector. \cr
+}
+
+}
+
+\subsection{Summarize inferences}{\tabular{ll}{
+   \strong{Method} \tab \strong{Description} \cr
+   \code{\link[=fit-method-summary]{$summary()}} \tab Run \code{\link[posterior:draws_summary]{posterior::summarise_draws()}}. \cr
+   \code{\link[=fit-method-cmdstan_summary]{$cmdstan_summary()}} \tab Run and print CmdStan's \code{bin/stansummary}. \cr
+}
+
+}
+
+\subsection{Save fitted model object and temporary files}{\tabular{ll}{
+   \strong{Method} \tab \strong{Description} \cr
+   \code{\link[=fit-method-save_object]{$save_object()}} \tab Save fitted model object to a file. \cr
+   \code{\link[=fit-method-save_output_files]{$save_output_files()}} \tab Save output CSV files to a specified location. \cr
+   \code{\link[=fit-method-save_data_file]{$save_data_file()}} \tab Save JSON data file to a specified location. \cr
+   \code{\link[=fit-method-save_latent_dynamics_files]{$save_latent_dynamics_files()}} \tab Save diagnostic CSV files to a specified location. \cr
+}
+
+}
+
+\subsection{Report run times, console output, return codes}{\tabular{ll}{
+   \strong{Method} \tab \strong{Description} \cr
+   \code{\link[=fit-method-time]{$time()}} \tab Report the total run time. \cr
+   \code{\link[=fit-method-output]{$output()}} \tab Pretty print the output that was printed to the console. \cr
+   \code{\link[=fit-method-return_codes]{$return_codes()}} \tab Return the return codes from the CmdStan runs. \cr
+}
+
+}
+}
+
+\seealso{
+The CmdStanR website
+(\href{https://mc-stan.org/cmdstanr/}{mc-stan.org/cmdstanr}) for online
+documentation and tutorials.
+
+The Stan and CmdStan documentation:
+\itemize{
+\item Stan documentation: \href{https://mc-stan.org/users/documentation/}{mc-stan.org/users/documentation}
+\item CmdStan User’s Guide: \href{https://mc-stan.org/docs/cmdstan-guide/}{mc-stan.org/docs/cmdstan-guide}
+}
+
+Other fitted model objects: 
+\code{\link{CmdStanDiagnose}},
+\code{\link{CmdStanGQ}},
+\code{\link{CmdStanLaplace}},
+\code{\link{CmdStanMCMC}},
+\code{\link{CmdStanMLE}},
+\code{\link{CmdStanVB}}
+}
+\concept{fitted model objects}
diff --git a/man/CmdStanVB.Rd b/man/CmdStanVB.Rd
index 9c8cd80b..323ec339 100644
--- a/man/CmdStanVB.Rd
+++ b/man/CmdStanVB.Rd
@@ -82,6 +82,7 @@ Other fitted model objects:
 \code{\link{CmdStanGQ}},
 \code{\link{CmdStanLaplace}},
 \code{\link{CmdStanMCMC}},
-\code{\link{CmdStanMLE}}
+\code{\link{CmdStanMLE}},
+\code{\link{CmdStanPathfinder}}
 }
 \concept{fitted model objects}
diff --git a/man/as_draws.CmdStanMCMC.Rd b/man/as_draws.CmdStanMCMC.Rd
index a63ae953..f28d5f50 100644
--- a/man/as_draws.CmdStanMCMC.Rd
+++ b/man/as_draws.CmdStanMCMC.Rd
@@ -7,6 +7,7 @@
 \alias{as_draws.CmdStanLaplace}
 \alias{as_draws.CmdStanVB}
 \alias{as_draws.CmdStanGQ}
+\alias{as_draws.CmdStanPathfinder}
 \title{Create a \code{draws} object from a CmdStanR fitted model object}
 \usage{
 \method{as_draws}{CmdStanMCMC}(x, ...)
@@ -18,6 +19,8 @@
 \method{as_draws}{CmdStanVB}(x, ...)
 
 \method{as_draws}{CmdStanGQ}(x, ...)
+
+\method{as_draws}{CmdStanPathfinder}(x, ...)
 }
 \arguments{
 \item{x}{A CmdStanR fitted model object.}
diff --git a/man/cmdstan_model.Rd b/man/cmdstan_model.Rd
index 6e9408dd..8ab7139a 100644
--- a/man/cmdstan_model.Rd
+++ b/man/cmdstan_model.Rd
@@ -104,13 +104,15 @@ fit_optim <- mod$optimize(data = my_data_file, jacobian = TRUE)
 fit_laplace <- mod$laplace(data = my_data_file, mode = fit_optim, draws = 2000)
 fit_laplace$summary()
 
-# Run 'variational' method to approximate the posterior (default is meanfield ADVI)
+# Run 'variational' method to use ADVI to approximate posterior
 fit_vb <- mod$variational(data = stan_data, seed = 123)
 fit_vb$summary()
-
-# Plot approximate posterior using bayesplot
 mcmc_hist(fit_vb$draws("theta"))
 
+# Run 'pathfinder' method, a new alternative to the variational method
+fit_pf <- mod$pathfinder(data = stan_data, seed = 123)
+fit_pf$summary()
+mcmc_hist(fit_pf$draws("theta"))
 
 # Specifying initial values as a function
 fit_mcmc_w_init_fun <- mod$sample(
diff --git a/man/cmdstanr-package.Rd b/man/cmdstanr-package.Rd
index 9f2feb18..aa695619 100644
--- a/man/cmdstanr-package.Rd
+++ b/man/cmdstanr-package.Rd
@@ -129,13 +129,15 @@ fit_optim <- mod$optimize(data = my_data_file, jacobian = TRUE)
 fit_laplace <- mod$laplace(data = my_data_file, mode = fit_optim, draws = 2000)
 fit_laplace$summary()
 
-# Run 'variational' method to approximate the posterior (default is meanfield ADVI)
+# Run 'variational' method to use ADVI to approximate posterior
 fit_vb <- mod$variational(data = stan_data, seed = 123)
 fit_vb$summary()
-
-# Plot approximate posterior using bayesplot
 mcmc_hist(fit_vb$draws("theta"))
 
+# Run 'pathfinder' method, a new alternative to the variational method
+fit_pf <- mod$pathfinder(data = stan_data, seed = 123)
+fit_pf$summary()
+mcmc_hist(fit_pf$draws("theta"))
 
 # Specifying initial values as a function
 fit_mcmc_w_init_fun <- mod$sample(
diff --git a/man/cmdstanr_example.Rd b/man/cmdstanr_example.Rd
index 3474811b..5efa2a42 100644
--- a/man/cmdstanr_example.Rd
+++ b/man/cmdstanr_example.Rd
@@ -7,7 +7,7 @@
 \usage{
 cmdstanr_example(
   example = c("logistic", "schools", "schools_ncp"),
-  method = c("sample", "optimize", "laplace", "variational", "diagnose"),
+  method = c("sample", "optimize", "laplace", "variational", "pathfinder", "diagnose"),
   ...,
   quiet = TRUE,
   force_recompile = getOption("cmdstanr_force_recompile", default = FALSE)
diff --git a/man/fit-method-lp.Rd b/man/fit-method-lp.Rd
index acda71b7..686c55d2 100644
--- a/man/fit-method-lp.Rd
+++ b/man/fit-method-lp.Rd
@@ -8,6 +8,8 @@
 \usage{
 lp()
 
+lp_approx()
+
 lp_approx()
 }
 \value{
diff --git a/man/model-method-check_syntax.Rd b/man/model-method-check_syntax.Rd
index 64e93c02..ab1fa8c3 100644
--- a/man/model-method-check_syntax.Rd
+++ b/man/model-method-check_syntax.Rd
@@ -85,6 +85,7 @@ Other CmdStanModel methods:
 \code{\link{model-method-generate-quantities}},
 \code{\link{model-method-laplace}},
 \code{\link{model-method-optimize}},
+\code{\link{model-method-pathfinder}},
 \code{\link{model-method-sample_mpi}},
 \code{\link{model-method-sample}},
 \code{\link{model-method-variables}},
diff --git a/man/model-method-compile.Rd b/man/model-method-compile.Rd
index 7e180574..52578854 100644
--- a/man/model-method-compile.Rd
+++ b/man/model-method-compile.Rd
@@ -156,6 +156,7 @@ Other CmdStanModel methods:
 \code{\link{model-method-generate-quantities}},
 \code{\link{model-method-laplace}},
 \code{\link{model-method-optimize}},
+\code{\link{model-method-pathfinder}},
 \code{\link{model-method-sample_mpi}},
 \code{\link{model-method-sample}},
 \code{\link{model-method-variables}},
diff --git a/man/model-method-diagnose.Rd b/man/model-method-diagnose.Rd
index 15e7e21c..413ad64b 100644
--- a/man/model-method-diagnose.Rd
+++ b/man/model-method-diagnose.Rd
@@ -127,6 +127,7 @@ Other CmdStanModel methods:
 \code{\link{model-method-generate-quantities}},
 \code{\link{model-method-laplace}},
 \code{\link{model-method-optimize}},
+\code{\link{model-method-pathfinder}},
 \code{\link{model-method-sample_mpi}},
 \code{\link{model-method-sample}},
 \code{\link{model-method-variables}},
diff --git a/man/model-method-expose_functions.Rd b/man/model-method-expose_functions.Rd
index 066c7bad..a62f7bb8 100644
--- a/man/model-method-expose_functions.Rd
+++ b/man/model-method-expose_functions.Rd
@@ -76,6 +76,7 @@ Other CmdStanModel methods:
 \code{\link{model-method-generate-quantities}},
 \code{\link{model-method-laplace}},
 \code{\link{model-method-optimize}},
+\code{\link{model-method-pathfinder}},
 \code{\link{model-method-sample_mpi}},
 \code{\link{model-method-sample}},
 \code{\link{model-method-variables}},
diff --git a/man/model-method-format.Rd b/man/model-method-format.Rd
index 294ad8cc..2aa34f18 100644
--- a/man/model-method-format.Rd
+++ b/man/model-method-format.Rd
@@ -105,6 +105,7 @@ Other CmdStanModel methods:
 \code{\link{model-method-generate-quantities}},
 \code{\link{model-method-laplace}},
 \code{\link{model-method-optimize}},
+\code{\link{model-method-pathfinder}},
 \code{\link{model-method-sample_mpi}},
 \code{\link{model-method-sample}},
 \code{\link{model-method-variables}},
diff --git a/man/model-method-generate-quantities.Rd b/man/model-method-generate-quantities.Rd
index adf2fd14..c9173c30 100644
--- a/man/model-method-generate-quantities.Rd
+++ b/man/model-method-generate-quantities.Rd
@@ -176,6 +176,7 @@ Other CmdStanModel methods:
 \code{\link{model-method-format}},
 \code{\link{model-method-laplace}},
 \code{\link{model-method-optimize}},
+\code{\link{model-method-pathfinder}},
 \code{\link{model-method-sample_mpi}},
 \code{\link{model-method-sample}},
 \code{\link{model-method-variables}},
diff --git a/man/model-method-laplace.Rd b/man/model-method-laplace.Rd
index 21e0e1b6..76f343be 100644
--- a/man/model-method-laplace.Rd
+++ b/man/model-method-laplace.Rd
@@ -151,13 +151,11 @@ deviation of the posterior distribution. If the mode is a maximum
 likelihood estimate (MLE), the sample provides an estimate of the standard
 error of the likelihood. Whether the mode is the MAP or MLE depends on
 the value of the \code{jacobian} argument when running optimization. See the
-\href{https://mc-stan.org/docs/cmdstan-guide/laplace-sampling.html}{Laplace Sampling}
-section of the CmdStan User's Guide for more details.
+\href{https://mc-stan.org/docs/cmdstan-guide/}{CmdStan User’s Guide}
+for more details.
 
 Any argument left as \code{NULL} will default to the default value used by the
-installed version of CmdStan. See the
-\href{https://mc-stan.org/docs/cmdstan-guide/}{CmdStan User’s Guide}
-for more details on the default arguments.
+installed version of CmdStan.
 }
 \examples{
 \dontrun{
@@ -199,6 +197,7 @@ Other CmdStanModel methods:
 \code{\link{model-method-format}},
 \code{\link{model-method-generate-quantities}},
 \code{\link{model-method-optimize}},
+\code{\link{model-method-pathfinder}},
 \code{\link{model-method-sample_mpi}},
 \code{\link{model-method-sample}},
 \code{\link{model-method-variables}},
diff --git a/man/model-method-optimize.Rd b/man/model-method-optimize.Rd
index 06535f14..6ce6bfe1 100644
--- a/man/model-method-optimize.Rd
+++ b/man/model-method-optimize.Rd
@@ -166,8 +166,8 @@ A \code{\link{CmdStanMLE}} object.
 The \verb{$optimize()} method of a \code{\link{CmdStanModel}} object runs
 Stan's optimizer to obtain a (penalized) maximum likelihood estimate or a
 maximum a posteriori estimate (if \code{jacobian=TRUE}). See the
-\href{https://mc-stan.org/docs/cmdstan-guide/maximum-likelihood-estimation.html}{Maximum Likelihood Estimation}
-section of the CmdStan User's Guide for more details.
+\href{https://mc-stan.org/docs/cmdstan-guide/index.html}{CmdStan User's Guide}
+for more details.
 
 Any argument left as \code{NULL} will default to the default value used by the
 installed version of CmdStan. See the \href{https://mc-stan.org/docs/cmdstan-guide/}{CmdStan User’s Guide} for more details on the
@@ -239,13 +239,15 @@ fit_optim <- mod$optimize(data = my_data_file, jacobian = TRUE)
 fit_laplace <- mod$laplace(data = my_data_file, mode = fit_optim, draws = 2000)
 fit_laplace$summary()
 
-# Run 'variational' method to approximate the posterior (default is meanfield ADVI)
+# Run 'variational' method to use ADVI to approximate posterior
 fit_vb <- mod$variational(data = stan_data, seed = 123)
 fit_vb$summary()
-
-# Plot approximate posterior using bayesplot
 mcmc_hist(fit_vb$draws("theta"))
 
+# Run 'pathfinder' method, a new alternative to the variational method
+fit_pf <- mod$pathfinder(data = stan_data, seed = 123)
+fit_pf$summary()
+mcmc_hist(fit_pf$draws("theta"))
 
 # Specifying initial values as a function
 fit_mcmc_w_init_fun <- mod$sample(
@@ -308,6 +310,7 @@ Other CmdStanModel methods:
 \code{\link{model-method-format}},
 \code{\link{model-method-generate-quantities}},
 \code{\link{model-method-laplace}},
+\code{\link{model-method-pathfinder}},
 \code{\link{model-method-sample_mpi}},
 \code{\link{model-method-sample}},
 \code{\link{model-method-variables}},
diff --git a/man/model-method-pathfinder.Rd b/man/model-method-pathfinder.Rd
new file mode 100644
index 00000000..3bc00d7a
--- /dev/null
+++ b/man/model-method-pathfinder.Rd
@@ -0,0 +1,327 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/model.R
+\name{model-method-pathfinder}
+\alias{model-method-pathfinder}
+\alias{pathfinder}
+\title{Run Stan's Pathfinder Variational Inference Algorithm}
+\usage{
+pathfinder(
+  data = NULL,
+  seed = NULL,
+  refresh = NULL,
+  init = NULL,
+  save_latent_dynamics = FALSE,
+  output_dir = NULL,
+  output_basename = NULL,
+  sig_figs = NULL,
+  opencl_ids = NULL,
+  num_threads = NULL,
+  init_alpha = NULL,
+  tol_obj = NULL,
+  tol_rel_obj = NULL,
+  tol_grad = NULL,
+  tol_rel_grad = NULL,
+  tol_param = NULL,
+  history_size = NULL,
+  single_path_draws = NULL,
+  draws = NULL,
+  num_paths = NULL,
+  max_lbfgs_iters = NULL,
+  num_elbo_draws = NULL,
+  save_single_paths = NULL
+)
+}
+\arguments{
+\item{data}{(multiple options) The data to use for the variables specified in
+the data block of the Stan program. One of the following:
+\itemize{
+\item A named list of \R objects with the names corresponding to variables
+declared in the data block of the Stan program. Internally this list is then
+written to JSON for CmdStan using \code{\link[=write_stan_json]{write_stan_json()}}. See
+\code{\link[=write_stan_json]{write_stan_json()}} for details on the conversions performed on \R objects
+before they are passed to Stan.
+\item A path to a data file compatible with CmdStan (JSON or \R dump). See the
+appendices in the CmdStan guide for details on using these formats.
+\item \code{NULL} or an empty list if the Stan program has no data block.
+}}
+
+\item{seed}{(positive integer(s)) A seed for the (P)RNG to pass to CmdStan.
+In the case of multi-chain sampling the single \code{seed} will automatically be
+augmented by the the run (chain) ID so that each chain uses a different
+seed. The exception is the transformed data block, which defaults to
+using same seed for all chains so that the same data is generated for all
+chains if RNG functions are used. The only time \code{seed} should be specified
+as a vector (one element per chain) is if RNG functions are used in
+transformed data and the goal is to generate \emph{different} data for each
+chain.}
+
+\item{refresh}{(non-negative integer) The number of iterations between
+printed screen updates. If \code{refresh = 0}, only error messages will be
+printed.}
+
+\item{init}{(multiple options) The initialization method to use for the
+variables declared in the parameters block of the Stan program. One of
+the following:
+\itemize{
+\item A real number \code{x>0}. This initializes \emph{all} parameters randomly between
+\verb{[-x,x]} on the \emph{unconstrained} parameter space.;
+\item The number \code{0}. This initializes \emph{all} parameters to \code{0};
+\item A character vector of paths (one per chain) to JSON or Rdump files
+containing initial values for all or some parameters. See
+\code{\link[=write_stan_json]{write_stan_json()}} to write \R objects to JSON files compatible with
+CmdStan.
+\item A list of lists containing initial values for all or some parameters. For
+MCMC the list should contain a sublist for each chain. For optimization and
+variational inference there should be just one sublist. The sublists should
+have named elements corresponding to the parameters for which you are
+specifying initial values. See \strong{Examples}.
+\item A function that returns a single list with names corresponding to the
+parameters for which you are specifying initial values. The function can
+take no arguments or a single argument \code{chain_id}. For MCMC, if the function
+has argument \code{chain_id} it will be supplied with the chain id (from 1 to
+number of chains) when called to generate the initial values. See
+\strong{Examples}.
+}}
+
+\item{save_latent_dynamics}{(logical) Should auxiliary diagnostic information
+about the latent dynamics be written to temporary diagnostic CSV files?
+This argument replaces CmdStan's \code{diagnostic_file} argument and the content
+written to CSV is controlled by the user's CmdStan installation and not
+CmdStanR (for some algorithms no content may be written). The default
+is \code{FALSE}, which is appropriate for almost every use case. To save the
+temporary files created when \code{save_latent_dynamics=TRUE} see the
+\code{\link[=fit-method-save_latent_dynamics_files]{$save_latent_dynamics_files()}}
+method.}
+
+\item{output_dir}{(string) A path to a directory where CmdStan should write
+its output CSV files. For interactive use this can typically be left at
+\code{NULL} (temporary directory) since CmdStanR makes the CmdStan output
+(posterior draws and diagnostics) available in \R via methods of the fitted
+model objects. The behavior of \code{output_dir} is as follows:
+\itemize{
+\item If \code{NULL} (the default), then the CSV files are written to a temporary
+directory and only saved permanently if the user calls one of the \verb{$save_*}
+methods of the fitted model object (e.g.,
+\code{\link[=fit-method-save_output_files]{$save_output_files()}}). These temporary
+files are removed when the fitted model object is
+\link[base:gc]{garbage collected} (manually or automatically).
+\item If a path, then the files are created in \code{output_dir} with names
+corresponding to the defaults used by \verb{$save_output_files()}.
+}}
+
+\item{output_basename}{(string) A string to use as a prefix for the names of
+the output CSV files of CmdStan. If \code{NULL} (the default), the basename of
+the output CSV files will be comprised from the model name, timestamp, and
+5 random characters.}
+
+\item{sig_figs}{(positive integer) The number of significant figures used
+when storing the output values. By default, CmdStan represent the output
+values with 6 significant figures. The upper limit for \code{sig_figs} is 18.
+Increasing this value will result in larger output CSV files and thus an
+increased usage of disk space.}
+
+\item{opencl_ids}{(integer vector of length 2) The platform and
+device IDs of the OpenCL device to use for fitting. The model must
+be compiled with \code{cpp_options = list(stan_opencl = TRUE)} for this
+argument to have an effect.}
+
+\item{num_threads}{(positive integer) If the model was
+\link[=model-method-compile]{compiled} with threading support, the number of
+threads to use in parallelized sections (e.g., for multi-path pathfinder
+as well as \code{reduce_sum}).}
+
+\item{init_alpha}{(positive real) The initial step size parameter.}
+
+\item{tol_obj}{(positive real) Convergence tolerance on changes in objective function value.}
+
+\item{tol_rel_obj}{(positive real) Convergence tolerance on relative changes in objective function value.}
+
+\item{tol_grad}{(positive real) Convergence tolerance on the norm of the gradient.}
+
+\item{tol_rel_grad}{(positive real) Convergence tolerance on the relative norm of the gradient.}
+
+\item{tol_param}{(positive real) Convergence tolerance on changes in parameter value.}
+
+\item{history_size}{(positive integer) The size of the history used when
+approximating the Hessian.}
+
+\item{single_path_draws}{(positive integer) Number of draws a single
+pathfinder should return. The number of draws PSIS sampling samples from
+will be equal to \code{single_path_draws * num_paths}.}
+
+\item{draws}{(positive integer) Number of draws to return after performing
+pareto smooted importance sampling (PSIS). This must be smaller than
+\code{single_path_draws * num_paths}.}
+
+\item{num_paths}{(positive integer) Number of single pathfinders to run.}
+
+\item{max_lbfgs_iters}{(positive integer) The maximum number of iterations
+for LBFGS.}
+
+\item{num_elbo_draws}{(positive integer) Number of draws to make when
+calculating the ELBO of the approximation at each iteration of LBFGS.}
+
+\item{save_single_paths}{(logical) Whether to save the results of single
+pathfinder runs in multi-pathfinder.}
+}
+\value{
+A \code{\link{CmdStanPathfinder}} object.
+}
+\description{
+The \verb{$pathfinder()} method of a \code{\link{CmdStanModel}} object runs
+Stan's Pathfinder algorithms. Pathfinder is a variational method for
+approximately sampling from differentiable log densities. Starting from a
+random initialization, Pathfinder locates normal approximations to the
+target density along a quasi-Newton optimization path, with local
+covariance estimated using the negative inverse Hessian estimates produced
+by the L-BFGS optimizer. Pathfinder returns draws from the Gaussian
+approximation with the lowest estimated Kullback-Leibler (KL) divergence to
+the true posterior. See the
+\href{https://mc-stan.org/docs/cmdstan-guide/}{CmdStan User’s Guide}
+for more details.
+
+Any argument left as \code{NULL} will default to the default value used by the
+installed version of CmdStan
+}
+\examples{
+\dontrun{
+library(cmdstanr)
+library(posterior)
+library(bayesplot)
+color_scheme_set("brightblue")
+
+# Set path to CmdStan
+# (Note: if you installed CmdStan via install_cmdstan() with default settings
+# then setting the path is unnecessary but the default below should still work.
+# Otherwise use the `path` argument to specify the location of your
+# CmdStan installation.)
+set_cmdstan_path(path = NULL)
+
+# Create a CmdStanModel object from a Stan program,
+# here using the example model that comes with CmdStan
+file <- file.path(cmdstan_path(), "examples/bernoulli/bernoulli.stan")
+mod <- cmdstan_model(file)
+mod$print()
+
+# Data as a named list (like RStan)
+stan_data <- list(N = 10, y = c(0,1,0,0,0,0,0,0,0,1))
+
+# Run MCMC using the 'sample' method
+fit_mcmc <- mod$sample(
+  data = stan_data,
+  seed = 123,
+  chains = 2,
+  parallel_chains = 2
+)
+
+# Use 'posterior' package for summaries
+fit_mcmc$summary()
+
+# Check sampling diagnostics
+fit_mcmc$diagnostic_summary()
+
+# Get posterior draws
+draws <- fit_mcmc$draws()
+print(draws)
+
+# Convert to data frame using posterior::as_draws_df
+as_draws_df(draws)
+
+# Plot posterior using bayesplot (ggplot2)
+mcmc_hist(fit_mcmc$draws("theta"))
+
+# For models fit using MCMC, if you like working with RStan's stanfit objects
+# then you can create one with rstan::read_stan_csv()
+# stanfit <- rstan::read_stan_csv(fit_mcmc$output_files())
+
+
+# Run 'optimize' method to get a point estimate (default is Stan's LBFGS algorithm)
+# and also demonstrate specifying data as a path to a file instead of a list
+my_data_file <- file.path(cmdstan_path(), "examples/bernoulli/bernoulli.data.json")
+fit_optim <- mod$optimize(data = my_data_file, seed = 123)
+fit_optim$summary()
+
+# Run 'optimize' again with 'jacobian=TRUE' and then draw from laplace approximation
+# to the posterior
+fit_optim <- mod$optimize(data = my_data_file, jacobian = TRUE)
+fit_laplace <- mod$laplace(data = my_data_file, mode = fit_optim, draws = 2000)
+fit_laplace$summary()
+
+# Run 'variational' method to use ADVI to approximate posterior
+fit_vb <- mod$variational(data = stan_data, seed = 123)
+fit_vb$summary()
+mcmc_hist(fit_vb$draws("theta"))
+
+# Run 'pathfinder' method, a new alternative to the variational method
+fit_pf <- mod$pathfinder(data = stan_data, seed = 123)
+fit_pf$summary()
+mcmc_hist(fit_pf$draws("theta"))
+
+# Specifying initial values as a function
+fit_mcmc_w_init_fun <- mod$sample(
+  data = stan_data,
+  seed = 123,
+  chains = 2,
+  refresh = 0,
+  init = function() list(theta = runif(1))
+)
+fit_mcmc_w_init_fun_2 <- mod$sample(
+  data = stan_data,
+  seed = 123,
+  chains = 2,
+  refresh = 0,
+  init = function(chain_id) {
+    # silly but demonstrates optional use of chain_id
+    list(theta = 1 / (chain_id + 1))
+  }
+)
+fit_mcmc_w_init_fun_2$init()
+
+# Specifying initial values as a list of lists
+fit_mcmc_w_init_list <- mod$sample(
+  data = stan_data,
+  seed = 123,
+  chains = 2,
+  refresh = 0,
+  init = list(
+    list(theta = 0.75), # chain 1
+    list(theta = 0.25)  # chain 2
+  )
+)
+fit_optim_w_init_list <- mod$optimize(
+  data = stan_data,
+  seed = 123,
+  init = list(
+    list(theta = 0.75)
+  )
+)
+fit_optim_w_init_list$init()
+}
+
+}
+\seealso{
+The CmdStanR website
+(\href{https://mc-stan.org/cmdstanr/}{mc-stan.org/cmdstanr}) for online
+documentation and tutorials.
+
+The Stan and CmdStan documentation:
+\itemize{
+\item Stan documentation: \href{https://mc-stan.org/users/documentation/}{mc-stan.org/users/documentation}
+\item CmdStan User’s Guide: \href{https://mc-stan.org/docs/cmdstan-guide/}{mc-stan.org/docs/cmdstan-guide}
+}
+
+Other CmdStanModel methods: 
+\code{\link{model-method-check_syntax}},
+\code{\link{model-method-compile}},
+\code{\link{model-method-diagnose}},
+\code{\link{model-method-expose_functions}},
+\code{\link{model-method-format}},
+\code{\link{model-method-generate-quantities}},
+\code{\link{model-method-laplace}},
+\code{\link{model-method-optimize}},
+\code{\link{model-method-sample_mpi}},
+\code{\link{model-method-sample}},
+\code{\link{model-method-variables}},
+\code{\link{model-method-variational}}
+}
+\concept{CmdStanModel methods}
diff --git a/man/model-method-sample.Rd b/man/model-method-sample.Rd
index 1c415f6a..ca42f1b8 100644
--- a/man/model-method-sample.Rd
+++ b/man/model-method-sample.Rd
@@ -348,13 +348,15 @@ fit_optim <- mod$optimize(data = my_data_file, jacobian = TRUE)
 fit_laplace <- mod$laplace(data = my_data_file, mode = fit_optim, draws = 2000)
 fit_laplace$summary()
 
-# Run 'variational' method to approximate the posterior (default is meanfield ADVI)
+# Run 'variational' method to use ADVI to approximate posterior
 fit_vb <- mod$variational(data = stan_data, seed = 123)
 fit_vb$summary()
-
-# Plot approximate posterior using bayesplot
 mcmc_hist(fit_vb$draws("theta"))
 
+# Run 'pathfinder' method, a new alternative to the variational method
+fit_pf <- mod$pathfinder(data = stan_data, seed = 123)
+fit_pf$summary()
+mcmc_hist(fit_pf$draws("theta"))
 
 # Specifying initial values as a function
 fit_mcmc_w_init_fun <- mod$sample(
@@ -418,6 +420,7 @@ Other CmdStanModel methods:
 \code{\link{model-method-generate-quantities}},
 \code{\link{model-method-laplace}},
 \code{\link{model-method-optimize}},
+\code{\link{model-method-pathfinder}},
 \code{\link{model-method-sample_mpi}},
 \code{\link{model-method-variables}},
 \code{\link{model-method-variational}}
diff --git a/man/model-method-sample_mpi.Rd b/man/model-method-sample_mpi.Rd
index 8c17a4de..1ffcbda5 100644
--- a/man/model-method-sample_mpi.Rd
+++ b/man/model-method-sample_mpi.Rd
@@ -316,6 +316,7 @@ Other CmdStanModel methods:
 \code{\link{model-method-generate-quantities}},
 \code{\link{model-method-laplace}},
 \code{\link{model-method-optimize}},
+\code{\link{model-method-pathfinder}},
 \code{\link{model-method-sample}},
 \code{\link{model-method-variables}},
 \code{\link{model-method-variational}}
diff --git a/man/model-method-variables.Rd b/man/model-method-variables.Rd
index b5b85f1c..dc80ed9a 100644
--- a/man/model-method-variables.Rd
+++ b/man/model-method-variables.Rd
@@ -45,6 +45,7 @@ Other CmdStanModel methods:
 \code{\link{model-method-generate-quantities}},
 \code{\link{model-method-laplace}},
 \code{\link{model-method-optimize}},
+\code{\link{model-method-pathfinder}},
 \code{\link{model-method-sample_mpi}},
 \code{\link{model-method-sample}},
 \code{\link{model-method-variational}}
diff --git a/man/model-method-variational.Rd b/man/model-method-variational.Rd
index b04f8b7c..019da21d 100644
--- a/man/model-method-variational.Rd
+++ b/man/model-method-variational.Rd
@@ -163,22 +163,17 @@ A \code{\link{CmdStanVB}} object.
 }
 \description{
 The \verb{$variational()} method of a \code{\link{CmdStanModel}} object runs
-Stan's variational Bayes (ADVI) algorithms.
-
-Any argument left as \code{NULL} will default to the default value used by the
-installed version of CmdStan. See the
+Stan's Automatic Differentiation Variational Inference (ADVI) algorithms.
+The approximation is a Gaussian in the unconstrained variable space. Stan
+implements two ADVI algorithms: the \code{algorithm="meanfield"} option uses a
+fully factorized Gaussian for the approximation; the \code{algorithm="fullrank"}
+option uses a Gaussian with a full-rank covariance matrix for the
+approximation. See the
 \href{https://mc-stan.org/docs/cmdstan-guide/}{CmdStan User’s Guide}
 for more details.
-}
-\details{
-CmdStan can fit a variational approximation to the posterior. The
-approximation is a Gaussian in the unconstrained variable space. Stan
-implements two variational algorithms. The \code{algorithm="meanfield"} option
-uses a fully factorized Gaussian for the approximation. The
-\code{algorithm="fullrank"} option uses a Gaussian with a full-rank covariance
-matrix for the approximation.
-
--- \href{https://github.com/stan-dev/cmdstan/releases/latest}{\emph{CmdStan Interface User's Guide}}
+
+Any argument left as \code{NULL} will default to the default value used by the
+installed version of CmdStan.
 }
 \examples{
 \dontrun{
@@ -244,13 +239,15 @@ fit_optim <- mod$optimize(data = my_data_file, jacobian = TRUE)
 fit_laplace <- mod$laplace(data = my_data_file, mode = fit_optim, draws = 2000)
 fit_laplace$summary()
 
-# Run 'variational' method to approximate the posterior (default is meanfield ADVI)
+# Run 'variational' method to use ADVI to approximate posterior
 fit_vb <- mod$variational(data = stan_data, seed = 123)
 fit_vb$summary()
-
-# Plot approximate posterior using bayesplot
 mcmc_hist(fit_vb$draws("theta"))
 
+# Run 'pathfinder' method, a new alternative to the variational method
+fit_pf <- mod$pathfinder(data = stan_data, seed = 123)
+fit_pf$summary()
+mcmc_hist(fit_pf$draws("theta"))
 
 # Specifying initial values as a function
 fit_mcmc_w_init_fun <- mod$sample(
@@ -314,6 +311,7 @@ Other CmdStanModel methods:
 \code{\link{model-method-generate-quantities}},
 \code{\link{model-method-laplace}},
 \code{\link{model-method-optimize}},
+\code{\link{model-method-pathfinder}},
 \code{\link{model-method-sample_mpi}},
 \code{\link{model-method-sample}},
 \code{\link{model-method-variables}}
diff --git a/tests/testthat/test-model-pathfinder.R b/tests/testthat/test-model-pathfinder.R
new file mode 100644
index 00000000..2f23cc9b
--- /dev/null
+++ b/tests/testthat/test-model-pathfinder.R
@@ -0,0 +1,144 @@
+context("model-pathfinder")
+
+set_cmdstan_path()
+stan_program <- testing_stan_file("bernoulli")
+mod <- testing_model("bernoulli")
+stan_program_fp <- testing_stan_file("bernoulli_fp")
+mod_fp <- testing_model("bernoulli_fp")
+
+# valid ways to supply data
+data_list <- testing_data("bernoulli")
+data_file_r <- test_path("resources", "data", "bernoulli.data.R")
+data_file_json <- test_path("resources", "data", "bernoulli.data.json")
+
+# these are all valid for sample()
+ok_arg_values <- list(
+  data = data_list,
+  output_dir = tempdir(),
+  refresh = 5,
+  init = 1.5,
+  seed = 12345,
+  init_alpha = 1,
+  tol_obj = 1e-12,
+  tol_rel_obj = 1e-12,
+  tol_grad = 1e-12,
+  tol_rel_grad = 1e-12,
+  tol_param = 1e-12,
+  history_size = 5,
+  num_elbo_draws = 10,
+  single_path_draws = 10,
+  draws = 100,
+  num_paths = 4,
+  max_lbfgs_iters = 100,
+  save_single_paths = FALSE)
+
+# using any one of these should cause sample() to error
+bad_arg_values <- list(
+  data = "NOT_A_FILE",
+  output_dir = "NOT_A_DIRECTORY",
+  refresh = -1,
+  init = "maybe :P",
+  seed = -80,
+  init_alpha = "cat.jpeg",
+  init_alpha = -3,
+  tol_obj = -1,
+  tol_rel_obj = -4,
+  tol_grad = -5,
+  tol_rel_grad = -9,
+  tol_param = -2,
+  history_size = -6,
+  num_elbo_draws = -8,
+  draws = "no thanks",
+  single_path_draws = "Just one plz",
+  num_paths = -1,
+  max_lbfgs_iters = "idk :/"
+)
+
+bad_arg_values_2 <- list(
+  data = "NOT_A_FILE",
+  output_dir = "NOT_A_DIRECTORY",
+  refresh = -1,
+  init = "maybe :P",
+  seed = -80,
+  init_alpha = -3,
+  tol_obj = -1,
+  tol_rel_obj = -4,
+  tol_grad = -5,
+  tol_rel_grad = -9,
+  tol_param = -2,
+  history_size = -6,
+  num_elbo_draws = -8,
+  draws = "no thanks",
+  single_path_draws = "nope",
+  num_paths = -1,
+  max_lbfgs_iters = "idk :/",
+  save_single_paths = "Mby"
+)
+
+bad_arg_values_3 <- list(
+  data = matrix(1:10),
+  output_dir = 8,
+  init = "maybe :P",
+  seed = -80,
+  init_alpha = -3,
+  tol_obj = -1,
+  tol_rel_obj = -4,
+  tol_grad = -5,
+  tol_rel_grad = -9,
+  tol_param = -2,
+  history_size = -6,
+  num_elbo_draws = -8,
+  draws = "no thanks",
+  single_path_draws = " ",
+  num_paths = "NO!",
+  max_lbfgs_iters = "idk :/"
+)
+expect_pathfinder_output <- function(object, num_chains = NULL) {
+  expect_output(object, regexp = "Finished in (.*) seconds.")
+}
+
+
+test_that("Pathfinder Runs", {
+  expect_pathfinder_output(fit <- mod$pathfinder(data=data_list, seed=1234, refresh = 0))
+  expect_is(fit, "CmdStanPathfinder")
+})
+
+test_that("pathfinder() method works with data files", {
+  expect_pathfinder_output(fit_r <- mod$pathfinder(data = data_file_r))
+  expect_is(fit_r, "CmdStanPathfinder")
+
+  expect_pathfinder_output(fit_json <- mod$pathfinder(data = data_file_json))
+  expect_is(fit_json, "CmdStanPathfinder")
+})
+
+test_that("pathfinder() method works with init file", {
+  init_list <- list(theta = 0.5)
+  init_file <- tempfile(
+    tmpdir = cmdstanr:::cmdstan_tempdir(),
+    pattern = "testing-inits-",
+    fileext = ".json"
+  )
+  write_stan_json(init_list, file = init_file)
+  expect_pathfinder_output(mod$pathfinder(data = data_file_r, init = init_file))
+})
+
+test_that("pathfinder() method runs when all arguments specified", {
+  expect_pathfinder_output(fit <- do.call(mod$pathfinder, ok_arg_values))
+  expect_is(fit, "CmdStanPathfinder")
+})
+
+test_that("pathfinder() method runs when the stan file is removed", {
+  stan_file_tmp <- tempfile(pattern = "tmp", fileext = ".stan")
+  file.copy(stan_program, stan_file_tmp)
+  mod_tmp <- cmdstan_model(stan_file_tmp)
+  file.remove(stan_file_tmp)
+  expect_pathfinder_output(
+    mod_tmp$pathfinder(data = data_list)
+  )
+})
+
+test_that("no error when checking estimates after failure", {
+  fit <- cmdstanr_example("schools", method = "pathfinder", seed = 123) # optim always fails for this
+  expect_silent(fit$summary()) # no error
+})
+