diff --git a/src/stan_math_backend/Cpp_str.ml b/src/stan_math_backend/Cpp_str.ml
new file mode 100644
index 0000000000..5d9061705b
--- /dev/null
+++ b/src/stan_math_backend/Cpp_str.ml
@@ -0,0 +1,59 @@
+open Bytes
+
+(* A version of stdlib Bytes.escaped but 
+  uses octal output, rather than default decimal,
+  for escapes like \123
+  This allows cpp to read them as literals
+*)
+let escaped_b s =
+  let n = ref 0 in
+  for i = 0 to length s - 1 do
+    n :=
+      !n
+      +
+      match unsafe_get s i with
+      | '\"' | '\\' | '\n' | '\t' | '\r' | '\b' -> 2
+      | ' ' .. '~' -> 1
+      | _ -> 4
+  done ;
+  if !n = length s then copy s
+  else
+    let s' = create !n in
+    n := 0 ;
+    for i = 0 to length s - 1 do
+      ( match unsafe_get s i with
+      | ('\"' | '\\') as c ->
+          unsafe_set s' !n '\\' ; incr n ; unsafe_set s' !n c
+      | '\n' -> unsafe_set s' !n '\\' ; incr n ; unsafe_set s' !n 'n'
+      | '\t' -> unsafe_set s' !n '\\' ; incr n ; unsafe_set s' !n 't'
+      | '\r' -> unsafe_set s' !n '\\' ; incr n ; unsafe_set s' !n 'r'
+      | '\b' -> unsafe_set s' !n '\\' ; incr n ; unsafe_set s' !n 'b'
+      | ' ' .. '~' as c -> unsafe_set s' !n c
+      | c ->
+          let a = Char.code c in
+          unsafe_set s' !n '\\' ;
+          incr n ;
+          (* changed *)
+          unsafe_set s' !n (Char.chr (48 + (a / 64))) ;
+          incr n ;
+          (* changed *)
+          unsafe_set s' !n (Char.chr (48 + (a / 8 mod 8))) ;
+          incr n ;
+          (* changed *)
+          unsafe_set s' !n (Char.chr (48 + (a mod 8))) ) ;
+      incr n
+    done ;
+    s'
+
+open String
+
+let escaped s : string =
+  let rec escape_if_needed (s : string) n i =
+    if i >= n then s
+    else
+      match unsafe_get s i with
+      | '\"' | '\\' | '\000' .. '\031' | '\127' .. '\255' ->
+          Bytes.to_string (escaped_b (Bytes.of_string s))
+      | _ -> escape_if_needed s n (i + 1)
+  in
+  escape_if_needed s (length s) 0
diff --git a/src/stan_math_backend/Cpp_str.mli b/src/stan_math_backend/Cpp_str.mli
new file mode 100644
index 0000000000..19b0e73085
--- /dev/null
+++ b/src/stan_math_backend/Cpp_str.mli
@@ -0,0 +1 @@
+val escaped : string -> string
diff --git a/src/stan_math_backend/Expression_gen.ml b/src/stan_math_backend/Expression_gen.ml
index 1ea2405122..734e6d0aa4 100644
--- a/src/stan_math_backend/Expression_gen.ml
+++ b/src/stan_math_backend/Expression_gen.ml
@@ -539,7 +539,7 @@ and pp_indexed_simple ppf (obj, idcs) =
 and pp_expr ppf Expr.Fixed.({pattern; meta} as e) =
   match pattern with
   | Var s -> pf ppf "%s" s
-  | Lit (Str, s) -> pf ppf "%S" s
+  | Lit (Str, s) -> pf ppf "\"%s\"" (Cpp_str.escaped s)
   | Lit (_, s) -> pf ppf "%s" s
   | FunApp
       ( StanLib (op, _, _)
diff --git a/test/integration/good/code-gen/cpp.expected b/test/integration/good/code-gen/cpp.expected
index 5d7b98ef64..05744db814 100644
--- a/test/integration/good/code-gen/cpp.expected
+++ b/test/integration/good/code-gen/cpp.expected
@@ -21265,6 +21265,326 @@ stan::math::profile_map& get_stan_profile_data() {
 
 
 
+  $ ../../../../../install/default/bin/stanc --print-cpp print_unicode.stan
+
+// Code generated by %%NAME%% %%VERSION%%
+#include <stan/model/model_header.hpp>
+namespace print_unicode_model_namespace {
+
+using stan::io::dump;
+using stan::model::assign;
+using stan::model::index_uni;
+using stan::model::index_max;
+using stan::model::index_min;
+using stan::model::index_min_max;
+using stan::model::index_multi;
+using stan::model::index_omni;
+using stan::model::model_base_crtp;
+using stan::model::rvalue;
+using namespace stan::math;
+
+
+stan::math::profile_map profiles__;
+static constexpr std::array<const char*, 3> locations_array__ = 
+{" (found before start of program)",
+ " (in 'print_unicode.stan', line 2, column 2 to column 24)",
+ " (in 'print_unicode.stan', line 3, column 2 to column 23)"};
+
+
+
+class print_unicode_model final : public model_base_crtp<print_unicode_model> {
+
+ private:
+   
+  
+ 
+ public:
+  ~print_unicode_model() { }
+  
+  inline std::string model_name() const final { return "print_unicode_model"; }
+
+  inline std::vector<std::string> model_compile_info() const noexcept {
+    return std::vector<std::string>{"stanc_version = %%NAME%%3 %%VERSION%%", "stancflags = --print-cpp"};
+  }
+  
+  
+  print_unicode_model(stan::io::var_context& context__,
+                      unsigned int random_seed__ = 0,
+                      std::ostream* pstream__ = nullptr) : model_base_crtp(0) {
+    int current_statement__ = 0;
+    using local_scalar_t__ = double ;
+    boost::ecuyer1988 base_rng__ = 
+        stan::services::util::create_rng(random_seed__, 0);
+    (void) base_rng__;  // suppress unused var warning
+    static constexpr const char* function__ = "print_unicode_model_namespace::print_unicode_model";
+    (void) function__;  // suppress unused var warning
+    local_scalar_t__ DUMMY_VAR__(std::numeric_limits<double>::quiet_NaN());
+    (void) DUMMY_VAR__;  // suppress unused var warning
+    try {
+      int pos__;
+      pos__ = 1;
+      current_statement__ = 1;
+      if (pstream__) {
+        stan_print(pstream__, "test: \320\211\360\237\230\203");
+        stan_print(pstream__, "\n");
+      }
+      current_statement__ = 2;
+      if (pstream__) {
+        stan_print(pstream__, "\316\273 \316\262 \316\266 \317\200");
+        stan_print(pstream__, "\n");
+      }
+    } catch (const std::exception& e) {
+      stan::lang::rethrow_located(e, locations_array__[current_statement__]);
+    }
+    num_params_r__ = 0U;
+    
+  }
+  
+  template <bool propto__, bool jacobian__ , typename VecR, typename VecI, 
+  stan::require_vector_like_t<VecR>* = nullptr, 
+  stan::require_vector_like_vt<std::is_integral, VecI>* = nullptr> 
+  inline stan::scalar_type_t<VecR> log_prob_impl(VecR& params_r__,
+                                                 VecI& params_i__,
+                                                 std::ostream* pstream__ = nullptr) const {
+    using T__ = stan::scalar_type_t<VecR>;
+    using local_scalar_t__ = T__;
+    T__ lp__(0.0);
+    stan::math::accumulator<T__> lp_accum__;
+    stan::io::deserializer<local_scalar_t__> in__(params_r__, params_i__);
+    int current_statement__ = 0;
+    local_scalar_t__ DUMMY_VAR__(std::numeric_limits<double>::quiet_NaN());
+    (void) DUMMY_VAR__;  // suppress unused var warning
+    static constexpr const char* function__ = "print_unicode_model_namespace::log_prob";
+    (void) function__;  // suppress unused var warning
+    
+    try {
+      
+    } catch (const std::exception& e) {
+      stan::lang::rethrow_located(e, locations_array__[current_statement__]);
+    }
+    lp_accum__.add(lp__);
+    return lp_accum__.sum();
+    } // log_prob_impl() 
+    
+  template <typename RNG, typename VecR, typename VecI, typename VecVar, 
+  stan::require_vector_like_vt<std::is_floating_point, VecR>* = nullptr, 
+  stan::require_vector_like_vt<std::is_integral, VecI>* = nullptr, 
+  stan::require_std_vector_vt<std::is_floating_point, VecVar>* = nullptr> 
+  inline void write_array_impl(RNG& base_rng__, VecR& params_r__,
+                               VecI& params_i__, VecVar& vars__,
+                               const bool emit_transformed_parameters__ = true,
+                               const bool emit_generated_quantities__ = true,
+                               std::ostream* pstream__ = nullptr) const {
+    using local_scalar_t__ = double;
+    stan::io::deserializer<local_scalar_t__> in__(params_r__, params_i__);
+    stan::io::serializer<local_scalar_t__> out__(vars__);
+    static constexpr bool propto__ = true;
+    (void) propto__;
+    double lp__ = 0.0;
+    (void) lp__;  // dummy to suppress unused var warning
+    int current_statement__ = 0; 
+    stan::math::accumulator<double> lp_accum__;
+    local_scalar_t__ DUMMY_VAR__(std::numeric_limits<double>::quiet_NaN());
+    constexpr bool jacobian__ = false;
+    (void) DUMMY_VAR__;  // suppress unused var warning
+    static constexpr const char* function__ = "print_unicode_model_namespace::write_array";
+    (void) function__;  // suppress unused var warning
+    
+    try {
+      if (logical_negation((primitive_value(emit_transformed_parameters__) ||
+            primitive_value(emit_generated_quantities__)))) {
+        return ;
+      } 
+      if (logical_negation(emit_generated_quantities__)) {
+        return ;
+      } 
+    } catch (const std::exception& e) {
+      stan::lang::rethrow_located(e, locations_array__[current_statement__]);
+    }
+    } // write_array_impl() 
+    
+  template <typename VecVar, typename VecI, 
+  stan::require_std_vector_t<VecVar>* = nullptr, 
+  stan::require_vector_like_vt<std::is_integral, VecI>* = nullptr> 
+  inline void transform_inits_impl(VecVar& params_r__, VecI& params_i__,
+                                   VecVar& vars__,
+                                   std::ostream* pstream__ = nullptr) const {
+    using local_scalar_t__ = double;
+    stan::io::deserializer<local_scalar_t__> in__(params_r__, params_i__);
+    stan::io::serializer<local_scalar_t__> out__(vars__);
+    int current_statement__ = 0;
+    local_scalar_t__ DUMMY_VAR__(std::numeric_limits<double>::quiet_NaN());
+    
+    try {
+      int pos__;
+      pos__ = 1;
+    } catch (const std::exception& e) {
+      stan::lang::rethrow_located(e, locations_array__[current_statement__]);
+    }
+    } // transform_inits_impl() 
+    
+  inline void get_param_names(std::vector<std::string>& names__) const {
+    
+    names__ = std::vector<std::string>{};
+    
+    } // get_param_names() 
+    
+  inline void get_dims(std::vector<std::vector<size_t>>& dimss__) const {
+    
+    dimss__ = std::vector<std::vector<size_t>>{};
+    
+    } // get_dims() 
+    
+  inline void constrained_param_names(
+                                      std::vector<std::string>& param_names__,
+                                      bool emit_transformed_parameters__ = true,
+                                      bool emit_generated_quantities__ = true) const
+    final {
+    
+    
+    if (emit_transformed_parameters__) {
+      
+    }
+    
+    if (emit_generated_quantities__) {
+      
+    }
+    
+    } // constrained_param_names() 
+    
+  inline void unconstrained_param_names(
+                                        std::vector<std::string>& param_names__,
+                                        bool emit_transformed_parameters__ = true,
+                                        bool emit_generated_quantities__ = true) const
+    final {
+    
+    
+    if (emit_transformed_parameters__) {
+      
+    }
+    
+    if (emit_generated_quantities__) {
+      
+    }
+    
+    } // unconstrained_param_names() 
+    
+  inline std::string get_constrained_sizedtypes() const {
+    
+    return std::string("[]");
+    
+    } // get_constrained_sizedtypes() 
+    
+  inline std::string get_unconstrained_sizedtypes() const {
+    
+    return std::string("[]");
+    
+    } // get_unconstrained_sizedtypes() 
+    
+  
+    // Begin method overload boilerplate
+    template <typename RNG>
+    inline void write_array(RNG& base_rng,
+                            Eigen::Matrix<double,Eigen::Dynamic,1>& params_r,
+                            Eigen::Matrix<double,Eigen::Dynamic,1>& vars,
+                            const bool emit_transformed_parameters = true,
+                            const bool emit_generated_quantities = true,
+                            std::ostream* pstream = nullptr) const {
+      const size_t num_params__ = 0;
+      const size_t num_transformed = 0;
+      const size_t num_gen_quantities = 0;
+      std::vector<double> vars_vec(num_params__
+       + (emit_transformed_parameters * num_transformed)
+       + (emit_generated_quantities * num_gen_quantities));
+      std::vector<int> params_i;
+      write_array_impl(base_rng, params_r, params_i, vars_vec,
+          emit_transformed_parameters, emit_generated_quantities, pstream);
+      vars = Eigen::Map<Eigen::Matrix<double,Eigen::Dynamic,1>>(
+        vars_vec.data(), vars_vec.size());
+    }
+
+    template <typename RNG>
+    inline void write_array(RNG& base_rng, std::vector<double>& params_r,
+                            std::vector<int>& params_i,
+                            std::vector<double>& vars,
+                            bool emit_transformed_parameters = true,
+                            bool emit_generated_quantities = true,
+                            std::ostream* pstream = nullptr) const {
+      const size_t num_params__ = 0;
+      const size_t num_transformed = 0;
+      const size_t num_gen_quantities = 0;
+      vars.resize(num_params__
+        + (emit_transformed_parameters * num_transformed)
+        + (emit_generated_quantities * num_gen_quantities));
+      write_array_impl(base_rng, params_r, params_i, vars, emit_transformed_parameters, emit_generated_quantities, pstream);
+    }
+
+    template <bool propto__, bool jacobian__, typename T_>
+    inline T_ log_prob(Eigen::Matrix<T_,Eigen::Dynamic,1>& params_r,
+                       std::ostream* pstream = nullptr) const {
+      Eigen::Matrix<int, -1, 1> params_i;
+      return log_prob_impl<propto__, jacobian__>(params_r, params_i, pstream);
+    }
+
+    template <bool propto__, bool jacobian__, typename T__>
+    inline T__ log_prob(std::vector<T__>& params_r,
+                        std::vector<int>& params_i,
+                        std::ostream* pstream = nullptr) const {
+      return log_prob_impl<propto__, jacobian__>(params_r, params_i, pstream);
+    }
+
+
+    inline void transform_inits(const stan::io::var_context& context,
+                         Eigen::Matrix<double, Eigen::Dynamic, 1>& params_r,
+                         std::ostream* pstream = nullptr) const final {
+      std::vector<double> params_r_vec(params_r.size());
+      std::vector<int> params_i;
+      transform_inits(context, params_i, params_r_vec, pstream);
+      params_r = Eigen::Map<Eigen::Matrix<double,Eigen::Dynamic,1>>(
+        params_r_vec.data(), params_r_vec.size());
+    }
+
+  inline void transform_inits(const stan::io::var_context& context,
+                              std::vector<int>& params_i,
+                              std::vector<double>& vars,
+                              std::ostream* pstream__ = nullptr) const {
+    const std::array<std::string, 0> names__ = std::array<std::string, 0>{};
+    
+     std::vector<double> params_r_flat__;
+     for (auto&& param_name__ : names__) {
+       const auto param_vec__ = context.vals_r(param_name__);
+       params_r_flat__.reserve(params_r_flat__.size() + param_vec__.size());
+       for (auto&& param_val__ : param_vec__) {
+         params_r_flat__.push_back(param_val__);
+       }
+     }
+    vars.resize(params_r_flat__.size());
+    transform_inits_impl(params_r_flat__, params_i, vars, pstream__);
+    } // transform_inits() 
+    
+};
+}
+using stan_model = print_unicode_model_namespace::print_unicode_model;
+
+#ifndef USING_R
+
+// Boilerplate
+stan::model::model_base& new_model(
+        stan::io::var_context& data_context,
+        unsigned int seed,
+        std::ostream* msg_stream) {
+  stan_model* m = new stan_model(data_context, seed, msg_stream);
+  return *m;
+}
+
+stan::math::profile_map& get_stan_profile_data() {
+  return print_unicode_model_namespace::profiles__;
+}
+
+#endif
+
+
+
   $ ../../../../../install/default/bin/stanc --print-cpp reduce_sum_m1.stan
 
 // Code generated by %%NAME%% %%VERSION%%
diff --git a/test/integration/good/code-gen/print_unicode.stan b/test/integration/good/code-gen/print_unicode.stan
new file mode 100644
index 0000000000..de3c5c42f5
--- /dev/null
+++ b/test/integration/good/code-gen/print_unicode.stan
@@ -0,0 +1,4 @@
+transformed data {
+  print("test: Љ😃");
+  print("λ β ζ π");
+}