`newton_solver` halts sampling

Description

I've been trying to implement a simple sampler in the generated quantities block using the inverse CDF method. First, I draw a uniform random variable on (0, 1) then use newton_solver to sample an x. The method seems to work well for a few iterations, but then sampling completely halts. @rok-cesnovar has narrowed this down to cases when the uniform random variable is close to 0.

Example

library(cmdstanr)
library(tidyverse)

model_code <- "
functions {
  real lawless_generalized_gamma_cdf(real x, real k, real mu, real sigma) {
    real w = (log(x) - mu) / sigma;
    return gamma_p(k, k * exp(1 / sqrt(k) * w));
  }
  vector system(vector y, real k, real mu, real sigma, real u) {
    return [lawless_generalized_gamma_cdf(y[1]| k, mu, sigma) - u]';
  }
}
data {
  real k;
  real mu;
  real sigma;
  real u;

  real scaling_step;
  real f_tol;
  int max_steps;
}
generated quantities {
  real solution = solve_newton_tol(system, [mu]', scaling_step, f_tol, max_steps, k, mu, sigma, u)[1];
}
"

model_file <- write_stan_file(model_code)
model <- cmdstan_model(model_file)


u <- 5e-5

fit <- model$sample(
  data = list(
    k = 0.96, mu = 4.84, sigma = 0.97, u = u,
    scaling_step = 10e-3, f_tol = 10e-6, max_steps = 1
  ),
  chains = 1,
  iter_sampling = 1,
  fixed_param = TRUE
)

I've also tried trying to solve on the log scale (using the log cdf and log(u)) and the problem persists. The sampling simply halts and does not continue, no error messages.

Current Version:

v4.7.0

Reproduce the supplied cmdstanr model, focusing on solve_newton_tol/newton_solver with u = 5e-5 and the shown solver settings. Trace the solver behavior for uniform values near zero and determine why sampling halts without an error; done means the example completes reliably and the regression is covered by an appropriate test.