diff --git a/src/Evolution/Systems/Cce/CMakeLists.txt b/src/Evolution/Systems/Cce/CMakeLists.txt
index 185cfcda27171..8ce2803830e39 100644
--- a/src/Evolution/Systems/Cce/CMakeLists.txt
+++ b/src/Evolution/Systems/Cce/CMakeLists.txt
@@ -12,6 +12,7 @@ spectre_target_sources(
   BoundaryData.cpp
   Equations.cpp
   GaugeTransformBoundaryData.cpp
+  KleinGordonSource.cpp
   LinearOperators.cpp
   LinearSolve.cpp
   NewmanPenrose.cpp
@@ -47,6 +48,7 @@ spectre_target_headers(
   SpecBoundaryData.hpp
   SwshDerivatives.hpp
   System.hpp
+  KleinGordonSource.hpp
   KleinGordonSystem.hpp
   Tags.hpp
   WorldtubeBufferUpdater.hpp
diff --git a/src/Evolution/Systems/Cce/Components/KleinGordonCharacteristicEvolution.hpp b/src/Evolution/Systems/Cce/Components/KleinGordonCharacteristicEvolution.hpp
index 168f40734475b..242789df33ed8 100644
--- a/src/Evolution/Systems/Cce/Components/KleinGordonCharacteristicEvolution.hpp
+++ b/src/Evolution/Systems/Cce/Components/KleinGordonCharacteristicEvolution.hpp
@@ -9,6 +9,7 @@
 #include "Evolution/Systems/Cce/Actions/InitializeKleinGordonVariables.hpp"
 #include "Evolution/Systems/Cce/Actions/PrecomputeKleinGordonSourceVariables.hpp"
 #include "Evolution/Systems/Cce/Components/CharacteristicEvolution.hpp"
+#include "Evolution/Systems/Cce/KleinGordonSource.hpp"
 #include "Evolution/Systems/Cce/KleinGordonSystem.hpp"
 #include "Parallel/GlobalCache.hpp"
 #include "Parallel/Local.hpp"
@@ -89,6 +90,10 @@ struct KleinGordonCharacteristicEvolution
                           Actions::CalculatePsi0AndDerivAtInnerBoundary,
                           tmpl::list<>>,
       Actions::PrecomputeKleinGordonSourceVariables,
+      tmpl::transform<
+          bondi_hypersurface_step_tags,
+          tmpl::bind<::Actions::MutateApply,
+                     tmpl::bind<ComputeKleinGordonSource, tmpl::_1>>>,
       tmpl::transform<bondi_hypersurface_step_tags,
                       tmpl::bind<hypersurface_computation, tmpl::_1>>,
       Actions::FilterSwshVolumeQuantity<Tags::BondiH>,
@@ -126,6 +131,10 @@ struct KleinGordonCharacteristicEvolution
                           Actions::CalculatePsi0AndDerivAtInnerBoundary,
                           tmpl::list<>>,
       Actions::PrecomputeKleinGordonSourceVariables,
+      tmpl::transform<
+          bondi_hypersurface_step_tags,
+          tmpl::bind<::Actions::MutateApply,
+                     tmpl::bind<ComputeKleinGordonSource, tmpl::_1>>>,
       tmpl::transform<bondi_hypersurface_step_tags,
                       tmpl::bind<hypersurface_computation, tmpl::_1>>,
       Actions::FilterSwshVolumeQuantity<Tags::BondiH>,
diff --git a/src/Evolution/Systems/Cce/KleinGordonSource.cpp b/src/Evolution/Systems/Cce/KleinGordonSource.cpp
new file mode 100644
index 0000000000000..7c157f3b1c255
--- /dev/null
+++ b/src/Evolution/Systems/Cce/KleinGordonSource.cpp
@@ -0,0 +1,57 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#include "Evolution/Systems/Cce/KleinGordonSource.hpp"
+
+namespace Cce {
+
+void ComputeKleinGordonSource<Tags::BondiBeta>::apply(
+    gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 0>>*> kg_source_beta,
+    const Scalar<SpinWeighted<ComplexDataVector, 0>>& dy_kg_psi,
+    const Scalar<SpinWeighted<ComplexDataVector, 0>>& one_minus_y) {
+  get(*kg_source_beta) = 2. * M_PI * get(one_minus_y) * square(get(dy_kg_psi));
+}
+
+void ComputeKleinGordonSource<Tags::BondiQ>::apply(
+    gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 1>>*> kg_source_q,
+    const Scalar<SpinWeighted<ComplexDataVector, 0>>& dy_kg_psi,
+    const Scalar<SpinWeighted<ComplexDataVector, 1>>& eth_kg_psi) {
+  get(*kg_source_q) = 16. * M_PI * get(eth_kg_psi) * get(dy_kg_psi);
+}
+
+void ComputeKleinGordonSource<Tags::BondiU>::apply(
+    gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 1>>*> kg_source_u,
+    const size_t l_max, const size_t number_of_radial_points) {
+  const size_t volume_size =
+      Spectral::Swsh::number_of_swsh_collocation_points(l_max) *
+      number_of_radial_points;
+  get(*kg_source_u) = SpinWeighted<ComplexDataVector, 1>{volume_size, 0.0};
+}
+
+void ComputeKleinGordonSource<Tags::BondiW>::apply(
+    gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 0>>*> kg_source_w,
+    const Scalar<SpinWeighted<ComplexDataVector, 0>>& exp_2_beta,
+    const Scalar<SpinWeighted<ComplexDataVector, 0>>& bondi_r,
+    const Scalar<SpinWeighted<ComplexDataVector, 0>>& bondi_k,
+    const Scalar<SpinWeighted<ComplexDataVector, 2>>& bondi_j,
+    const Scalar<SpinWeighted<ComplexDataVector, 1>>& eth_kg_psi) {
+  SpinWeighted<ComplexDataVector, 0> complex_part;
+  SpinWeighted<ComplexDataVector, 0> real_part;
+
+  complex_part.data() =
+      2.0 * real(get(bondi_j).data() * square(conj(get(eth_kg_psi))).data());
+  real_part.data() =
+      -2.0 * get(bondi_k).data() * square(abs(get(eth_kg_psi).data()));
+  get(*kg_source_w) =
+      M_PI * get(exp_2_beta) / get(bondi_r) * (complex_part + real_part);
+}
+
+void ComputeKleinGordonSource<Tags::BondiH>::apply(
+    gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 2>>*> kg_source_h,
+    const Scalar<SpinWeighted<ComplexDataVector, 0>>& exp_2_beta,
+    const Scalar<SpinWeighted<ComplexDataVector, 0>>& bondi_r,
+    const Scalar<SpinWeighted<ComplexDataVector, 1>>& eth_kg_psi) {
+  get(*kg_source_h) =
+      2 * M_PI * get(exp_2_beta) / get(bondi_r) * square(get(eth_kg_psi));
+}
+}  // namespace Cce
diff --git a/src/Evolution/Systems/Cce/KleinGordonSource.hpp b/src/Evolution/Systems/Cce/KleinGordonSource.hpp
new file mode 100644
index 0000000000000..cf63916b632b4
--- /dev/null
+++ b/src/Evolution/Systems/Cce/KleinGordonSource.hpp
@@ -0,0 +1,139 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#pragma once
+
+#include "Evolution/Systems/Cce/OptionTags.hpp"
+#include "Evolution/Systems/Cce/Tags.hpp"
+#include "Utilities/Gsl.hpp"
+#include "Utilities/TMPL.hpp"
+
+namespace Cce {
+
+/*!
+ * \brief Computes `Tags::KleinGordonSource<Tag>` for the tags evolved by
+ * Klein-Gordon Cce.
+ *
+ * \details In scalar-tensor theory, the Cce hypersurface equations get
+ * additional source terms contributed by the stress-energy tensor of the scalar
+ * field. The tag `Tags::KleinGordonSource<Tag>` stores the corresponding volume
+ * data.
+ */
+template <typename Tag>
+struct ComputeKleinGordonSource;
+
+/*!
+ * \brief Computes the Klein-Gordon source of the Bondi \f$\beta\f$
+ *
+ * \details The source reads:
+ *
+ * \f{align*}{
+ * 2 \pi (1-y) (\partial_y\psi)^2,
+ * \f}
+ * where \f$\psi\f$ is the Klein-Gordon (scalar) field.
+ */
+template <>
+struct ComputeKleinGordonSource<Tags::BondiBeta> {
+  using return_tags = tmpl::list<Tags::KleinGordonSource<Tags::BondiBeta>>;
+  using argument_tags =
+      tmpl::list<Tags::Dy<Tags::KleinGordonPsi>, Tags::OneMinusY>;
+  static void apply(
+      gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 0>>*> kg_source_beta,
+      const Scalar<SpinWeighted<ComplexDataVector, 0>>& dy_kg_psi,
+      const Scalar<SpinWeighted<ComplexDataVector, 0>>& one_minus_y);
+};
+
+/*!
+ * \brief Computes the Klein-Gordon source of the Bondi \f$Q\f$
+ *
+ * \details Following the nomenclature of \cite Moxon2020gha and their Eq. (49),
+ * the scalar field contributes only to the regular part of the source term
+ * \f$S_2^R\f$. The expression reads:
+ *
+ * \f{align*}{
+ * 16 \pi \eth\psi \partial_y\psi,
+ * \f}
+ * where \f$\psi\f$ is the Klein-Gordon (scalar) field.
+ */
+template <>
+struct ComputeKleinGordonSource<Tags::BondiQ> {
+  using return_tags = tmpl::list<Tags::KleinGordonSource<Tags::BondiQ>>;
+  using argument_tags =
+      tmpl::list<Tags::Dy<Tags::KleinGordonPsi>,
+                 Spectral::Swsh::Tags::Derivative<Tags::KleinGordonPsi,
+                                                  Spectral::Swsh::Tags::Eth>>;
+  static void apply(
+      gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 1>>*> kg_source_q,
+      const Scalar<SpinWeighted<ComplexDataVector, 0>>& dy_kg_psi,
+      const Scalar<SpinWeighted<ComplexDataVector, 1>>& eth_kg_psi);
+};
+
+/*!
+ * \brief Computes the Klein-Gordon source of the Bondi \f$U\f$
+ *
+ * \details The source vanishes.
+ */
+template <>
+struct ComputeKleinGordonSource<Tags::BondiU> {
+  using return_tags = tmpl::list<Tags::KleinGordonSource<Tags::BondiU>>;
+  using argument_tags = tmpl::list<Tags::LMax, Tags::NumberOfRadialPoints>;
+  static void apply(
+      gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 1>>*> kg_source_u,
+      size_t l_max, size_t number_of_radial_points);
+};
+
+/*!
+ * \brief Computes the Klein-Gordon source of the Bondi \f$W\f$
+ *
+ * \details Following the nomenclature of \cite Moxon2020gha and their Eq. (49),
+ * the scalar field contributes only to the regular part of the source term
+ * \f$S_2^R\f$. The expression reads:
+ *
+ * \f{align*}{
+ * \frac{\pi e^{2\beta}}{R} \left[J(\bar{\eth}\psi)^2 + \bar{J}(\eth
+ * \psi)^2-2K \eth\psi\bar{\eth}\psi\right],
+ * \f}
+ * where \f$\psi\f$ is the Klein-Gordon (scalar) field.
+ */
+template <>
+struct ComputeKleinGordonSource<Tags::BondiW> {
+  using return_tags = tmpl::list<Tags::KleinGordonSource<Tags::BondiW>>;
+  using argument_tags =
+      tmpl::list<Tags::Exp2Beta, Tags::BondiR, Tags::BondiK, Tags::BondiJ,
+                 Spectral::Swsh::Tags::Derivative<Tags::KleinGordonPsi,
+                                                  Spectral::Swsh::Tags::Eth>>;
+  static void apply(
+      gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 0>>*> kg_source_w,
+      const Scalar<SpinWeighted<ComplexDataVector, 0>>& exp_2_beta,
+      const Scalar<SpinWeighted<ComplexDataVector, 0>>& bondi_r,
+      const Scalar<SpinWeighted<ComplexDataVector, 0>>& bondi_k,
+      const Scalar<SpinWeighted<ComplexDataVector, 2>>& bondi_j,
+      const Scalar<SpinWeighted<ComplexDataVector, 1>>& eth_kg_psi);
+};
+
+/*!
+ * \brief Computes the Klein-Gordon source of the Bondi \f$H\f$
+ *
+ * \details Following the nomenclature of \cite Moxon2020gha and their Eq. (50),
+ * the scalar field contributes only to the regular part of the source term
+ * \f$S_3^R\f$. The expression reads:
+ *
+ * \f{align*}{
+ * 2 \pi \frac{e^{2\beta}}{R}(\eth\psi)^2,
+ * \f}
+ * where \f$\psi\f$ is the Klein-Gordon (scalar) field.
+ */
+template <>
+struct ComputeKleinGordonSource<Tags::BondiH> {
+  using return_tags = tmpl::list<Tags::KleinGordonSource<Tags::BondiH>>;
+  using argument_tags =
+      tmpl::list<Tags::Exp2Beta, Tags::BondiR,
+                 Spectral::Swsh::Tags::Derivative<Tags::KleinGordonPsi,
+                                                  Spectral::Swsh::Tags::Eth>>;
+  static void apply(
+      gsl::not_null<Scalar<SpinWeighted<ComplexDataVector, 2>>*> kg_source_h,
+      const Scalar<SpinWeighted<ComplexDataVector, 0>>& exp_2_beta,
+      const Scalar<SpinWeighted<ComplexDataVector, 0>>& bondi_r,
+      const Scalar<SpinWeighted<ComplexDataVector, 1>>& eth_kg_psi);
+};
+}  // namespace Cce
diff --git a/tests/Unit/Evolution/Systems/Cce/Actions/CMakeLists.txt b/tests/Unit/Evolution/Systems/Cce/Actions/CMakeLists.txt
index bfddede97214c..c6bfa668465f1 100644
--- a/tests/Unit/Evolution/Systems/Cce/Actions/CMakeLists.txt
+++ b/tests/Unit/Evolution/Systems/Cce/Actions/CMakeLists.txt
@@ -17,6 +17,7 @@ set(LIBRARY_SOURCES
   Test_InitializeKleinGordonFirstHypersurface.cpp
   Test_InitializeKleinGordonWorldtubeBoundary.cpp
   Test_InitializeWorldtubeBoundary.cpp
+  Test_KleinGordonCceCalculations.cpp
   Test_KleinGordonH5BoundaryCommunication.cpp
   Test_Psi0Matching.cpp
   Test_RequestBoundaryData.cpp
diff --git a/tests/Unit/Evolution/Systems/Cce/Actions/Test_KleinGordonCceCalculations.cpp b/tests/Unit/Evolution/Systems/Cce/Actions/Test_KleinGordonCceCalculations.cpp
new file mode 100644
index 0000000000000..1e4ab27572602
--- /dev/null
+++ b/tests/Unit/Evolution/Systems/Cce/Actions/Test_KleinGordonCceCalculations.cpp
@@ -0,0 +1,159 @@
+// Distributed under the MIT License.
+// See LICENSE.txt for details.
+
+#include "Framework/TestingFramework.hpp"
+
+#include "Evolution/Systems/Cce/Actions/InitializeKleinGordonFirstHypersurface.hpp"
+#include "Evolution/Systems/Cce/KleinGordonSource.hpp"
+#include "Evolution/Systems/Cce/PrecomputeCceDependencies.hpp"
+#include "Framework/ActionTesting.hpp"
+#include "Framework/TestHelpers.hpp"
+#include "Helpers/DataStructures/MakeWithRandomValues.hpp"
+#include "Helpers/NumericalAlgorithms/Spectral/SwshTestHelpers.hpp"
+#include "NumericalAlgorithms/Spectral/SwshFiltering.hpp"
+#include "ParallelAlgorithms/Actions/MutateApply.hpp"
+#include "Utilities/Gsl.hpp"
+
+namespace Cce {
+namespace {
+
+using real_volume_tags_to_generate = tmpl::list<Tags::OneMinusY>;
+
+using swsh_volume_tags_to_generate =
+    tmpl::list<Cce::Tags::Exp2Beta, Cce::Tags::BondiR, Cce::Tags::BondiK,
+               Cce::Tags::BondiJ, Tags::Dy<Tags::KleinGordonPsi>,
+               Spectral::Swsh::Tags::Derivative<Cce::Tags::KleinGordonPsi,
+                                                Spectral::Swsh::Tags::Eth>>;
+
+using swsh_volume_tags_to_compute =
+    tmpl::list<Tags::KleinGordonSource<Tags::BondiBeta>,
+               Tags::KleinGordonSource<Tags::BondiQ>,
+               Tags::KleinGordonSource<Tags::BondiU>,
+               Tags::KleinGordonSource<Tags::BondiW>,
+               Tags::KleinGordonSource<Tags::BondiH>>;
+
+template <typename Metavariables>
+struct mock_kg_characteristic_evolution {
+  using simple_tags = db::AddSimpleTags<::Tags::Variables<
+      tmpl::append<real_volume_tags_to_generate, swsh_volume_tags_to_generate,
+                   swsh_volume_tags_to_compute>>>;
+
+  using metavariables = Metavariables;
+  using chare_type = ActionTesting::MockArrayChare;
+  using array_index = size_t;
+
+  using phase_dependent_action_list = tmpl::list<
+      Parallel::PhaseActions<
+          Parallel::Phase::Initialization,
+          tmpl::list<ActionTesting::InitializeDataBox<simple_tags>>>,
+      Parallel::PhaseActions<
+          Parallel::Phase::Evolve,
+          tmpl::list<tmpl::transform<
+              bondi_hypersurface_step_tags,
+              tmpl::bind<::Actions::MutateApply,
+                         tmpl::bind<ComputeKleinGordonSource, tmpl::_1>>>>>>;
+
+  using const_global_cache_tags =
+      tmpl::list<Tags::LMax, Tags::NumberOfRadialPoints>;
+};
+
+struct metavariables {
+  using component_list =
+      tmpl::list<mock_kg_characteristic_evolution<metavariables>>;
+};
+
+// This unit test is to validate the automatic calling chain of Klein-Gordon Cce
+// calculations, including the mutators `ComputeKleinGordonSource`. The test
+// involves
+// (a) Fills a bunch of variables with random numbers (filtered so that there is
+// no aliasing in highest modes).
+// (b) Puts those variables in two places: the MockRuntimeSystem runner and a
+// databox called expected_box.
+// (c) Calls next_action a few times (which fills the results in the databox of
+// MockRuntimeSystem runner).
+// (d) Calls all the individual calculations by hand, filling the results in
+// expected_box.
+// (e) Compares MockRuntimeSystem's databox vs expected_box.
+template <typename Generator>
+void test_klein_gordon_cce_source(const gsl::not_null<Generator*> gen) {
+  UniformCustomDistribution<size_t> sdist{7, 10};
+  const size_t l_max = sdist(*gen);
+  const size_t number_of_radial_points = sdist(*gen);
+  UniformCustomDistribution<double> coefficient_distribution{-2.0, 2.0};
+
+  Variables<
+      tmpl::append<real_volume_tags_to_generate, swsh_volume_tags_to_generate,
+                   swsh_volume_tags_to_compute>>
+      component_volume_variables{
+          Spectral::Swsh::number_of_swsh_collocation_points(l_max) *
+          number_of_radial_points};
+
+  tmpl::for_each<swsh_volume_tags_to_generate>(
+      [&component_volume_variables, &gen, &coefficient_distribution,
+       &number_of_radial_points, &l_max](auto tag_v) {
+        using tag = typename decltype(tag_v)::type;
+        SpinWeighted<ComplexModalVector, tag::type::type::spin> generated_modes{
+            number_of_radial_points *
+            Spectral::Swsh::size_of_libsharp_coefficient_vector(l_max)};
+        Spectral::Swsh::TestHelpers::generate_swsh_modes<tag::type::type::spin>(
+            make_not_null(&generated_modes.data()), gen,
+            make_not_null(&coefficient_distribution), number_of_radial_points,
+            l_max);
+        get(get<tag>(component_volume_variables)) =
+            Spectral::Swsh::inverse_swsh_transform(
+                l_max, number_of_radial_points, generated_modes);
+        // aggressive filter to make the uniformly generated random modes
+        // somewhat reasonable
+        Spectral::Swsh::filter_swsh_volume_quantity(
+            make_not_null(&get(get<tag>(component_volume_variables))), l_max,
+            l_max / 2, 32.0, 8);
+      });
+
+  PrecomputeCceDependencies<Tags::EvolutionGaugeBoundaryValue,
+                            Tags::OneMinusY>::
+      apply(make_not_null(&get<Tags::OneMinusY>(component_volume_variables)),
+            l_max, number_of_radial_points);
+
+  using component = mock_kg_characteristic_evolution<metavariables>;
+
+  // tests start here
+  ActionTesting::MockRuntimeSystem<metavariables> runner{
+      tuples::tagged_tuple_from_typelist<
+          Parallel::get_const_global_cache_tags<metavariables>>{
+          l_max, number_of_radial_points}};
+  ActionTesting::emplace_component_and_initialize<component>(
+      &runner, 0, {component_volume_variables});
+  auto expected_box = db::create<
+      tmpl::append<component::simple_tags,
+                   db::AddSimpleTags<Tags::LMax, Tags::NumberOfRadialPoints>>>(
+      component_volume_variables, l_max, number_of_radial_points);
+
+  runner.set_phase(Parallel::Phase::Evolve);
+
+  for (int i = 0; i < 5; i++) {
+    ActionTesting::next_action<component>(make_not_null(&runner), 0);
+  }
+
+  // tests for `ComputeKleinGordonSource`
+  tmpl::for_each<bondi_hypersurface_step_tags>([&expected_box,
+                                                &runner](auto tag_v) {
+    using tag = typename decltype(tag_v)::type;
+    db::mutate_apply<ComputeKleinGordonSource<tag>>(
+        make_not_null(&expected_box));
+
+    auto computed_result =
+        ActionTesting::get_databox_tag<component, Tags::KleinGordonSource<tag>>(
+            runner, 0);
+
+    auto expected_result = db::get<Tags::KleinGordonSource<tag>>(expected_box);
+    CHECK(computed_result == expected_result);
+  });
+}
+}  // namespace
+
+SPECTRE_TEST_CASE("Unit.Evolution.Systems.Cce.Actions.KGCceCalculations",
+                  "[Unit][Cce]") {
+  MAKE_GENERATOR(gen);
+  test_klein_gordon_cce_source(make_not_null(&gen));
+}
+}  // namespace Cce