diff --git a/docs/src/bath_boson/bosonic_bath_intro.md b/docs/src/bath_boson/bosonic_bath_intro.md
index 61f32fdb..a51373ff 100644
--- a/docs/src/bath_boson/bosonic_bath_intro.md
+++ b/docs/src/bath_boson/bosonic_bath_intro.md
@@ -56,9 +56,9 @@ BosonBath object with 4 exponential-expansion terms
 ```
 
 ## Calculate the correlation function
-To check whether the exponential terms in the [`BosonBath`](@ref) is correct or not, one can call [`C(bath::BosonBath, tlist::AbstractVector)`](@ref) to calculate the correlation function ``C(t)``, where ``t=t_1-t_2``:
+To check whether the exponential terms in the [`BosonBath`](@ref) is correct or not, one can call [`correlation_function`](@ref) to calculate the correlation function ``C(t)``, where ``t=t_1-t_2``:
 ```julia
-c_list = C(bath, tlist)
+c_list = correlation_function(bath, tlist)
 ```
 Here, `c_list` is a list which contains the value of ``C(t)`` corresponds to the given time series `tlist`.
 
diff --git a/docs/src/bath_boson_RWA/bosonic_bath_RWA_intro.md b/docs/src/bath_boson_RWA/bosonic_bath_RWA_intro.md
index d23b8406..27dc993a 100644
--- a/docs/src/bath_boson_RWA/bosonic_bath_RWA_intro.md
+++ b/docs/src/bath_boson_RWA/bosonic_bath_RWA_intro.md
@@ -41,9 +41,9 @@ BosonBath object with 4 exponential-expansion terms
 Note that [`BosonBath`](@ref) under RWA always have even number of exponential terms (half for ``C^{\nu=+}`` and half for ``C^{\nu=-}``)
 
 ## Calculate the correlation function
-To check whether the exponential terms in the [`FermionBath`](@ref) is correct or not, one can call [`C(bath::BosonBath, tlist::AbstractVector)`](@ref) to calculate the correlation function ``C(t)``, where ``t=t_1-t_2``:
+To check whether the exponential terms in the [`FermionBath`](@ref) is correct or not, one can call [`correlation_function`](@ref) to calculate the correlation function ``C(t)``, where ``t=t_1-t_2``:
 ```julia
-cp_list, cm_list = C(bath, tlist)
+cp_list, cm_list = correlation_function(bath, tlist)
 ```
 Here, `cp_list` and `cm_list` are the lists which contain the value of ``C^{\nu=+}(t)`` and ``C^{\nu=-}(t)`` correspond to the given time series `tlist`, respectively.
 
diff --git a/docs/src/bath_fermion/fermionic_bath_intro.md b/docs/src/bath_fermion/fermionic_bath_intro.md
index 15bc536c..e85b87e7 100644
--- a/docs/src/bath_fermion/fermionic_bath_intro.md
+++ b/docs/src/bath_fermion/fermionic_bath_intro.md
@@ -44,9 +44,9 @@ FermionBath object with 4 exponential-expansion terms
 Note that [`FermionBath`](@ref) always have even number of exponential terms (half for ``C^{\nu=+}`` and half for ``C^{\nu=-}``)
 
 ## Calculate the correlation function
-To check whether the exponential terms in the [`FermionBath`](@ref) is correct or not, one can call [`C(bath::FermionBath, tlist::AbstractVector)`](@ref) to calculate the correlation function ``C(t)``, where ``t=t_1-t_2``:
+To check whether the exponential terms in the [`FermionBath`](@ref) is correct or not, one can call [`correlation_function`](@ref) to calculate the correlation function ``C(t)``, where ``t=t_1-t_2``:
 ```julia
-cp_list, cm_list = C(bath, tlist)
+cp_list, cm_list = correlation_function(bath, tlist)
 ```
 Here, `cp_list` and `cm_list` are the lists which contain the value of ``C^{\nu=+}(t)`` and ``C^{\nu=-}(t)`` correspond to the given time series `tlist`, respectively.
 
diff --git a/docs/src/libraryAPI.md b/docs/src/libraryAPI.md
index ff276ef8..d0c8d7a7 100644
--- a/docs/src/libraryAPI.md
+++ b/docs/src/libraryAPI.md
@@ -15,8 +15,7 @@ Pages = ["libraryAPI.md"]
 ## [Bath Module](@id lib-Bath)
 
 ```@docs
-C(bath::BosonBath, tlist::AbstractVector)
-C(bath::FermionBath, tlist::AbstractVector)
+correlation_function
 Exponent
 BosonBath
 BosonBath(op::QuantumObject, η::Vector{Ti}, γ::Vector{Tj}, δ::Number=0.0; combine::Bool=true) where {Ti, Tj <: Number}
diff --git a/examples/cavityQED.jl b/examples/cavityQED.jl
index 768b0ac0..c01e3b2e 100644
--- a/examples/cavityQED.jl
+++ b/examples/cavityQED.jl
@@ -70,13 +70,13 @@ kT = 0.025
 N = 20
 Bath = Boson_DrudeLorentz_Pade(a + a', Γ, W, kT, N)
 
-# Before incorporating the correlation function into the HEOMLS matrix, it is essential to verify if the total number of exponentials for the reservoir sufficiently describes the practical situation.
+# Before incorporating the correlation function into the HEOMLS matrix, it is essential to verify (by using [`correlation_function`](@ref)) if the total number of exponentials for the reservoir sufficiently describes the practical situation.
 
 tlist_test = 0:0.1:10;
 
 Bath_test = Boson_DrudeLorentz_Pade(a + a', Γ, W, kT, 1000);
-Ct = C(Bath, tlist_test);
-Ct2 = C(Bath_test, tlist_test);
+Ct = correlation_function(Bath, tlist_test);
+Ct2 = correlation_function(Bath_test, tlist_test);
 
 Plots.plot(tlist_test, real(Ct), label = "N=20 (real part )", linestyle = :dash, linewidth = 3)
 Plots.plot!(tlist_test, real(Ct2), label = "N=1000 (real part)", linestyle = :solid, linewidth = 3)
diff --git a/src/bath/BathBase.jl b/src/bath/BathBase.jl
index 875bbe90..f0f34416 100644
--- a/src/bath/BathBase.jl
+++ b/src/bath/BathBase.jl
@@ -1,4 +1,4 @@
-export AbstractBath, Exponent, C
+export AbstractBath, Exponent, correlation_function
 
 abstract type AbstractBath end
 
@@ -145,3 +145,5 @@ function _check_gamma_absorb_and_emit(γ_absorb, γ_emit)
         error("The length of \'γ_absorb\' and \'γ_emit\' should be the same.")
     end
 end
+
+correlation_function(bath::AbstractBath, t::Real) = correlation_function(bath, [t])
diff --git a/src/bath/BosonBath.jl b/src/bath/BosonBath.jl
index 2b89c456..fd0f0d6a 100644
--- a/src/bath/BosonBath.jl
+++ b/src/bath/BosonBath.jl
@@ -496,7 +496,7 @@ function bosonEmit(
 end
 
 @doc raw"""
-    C(bath, tlist)
+    correlation_function(bath, tlist)
 Calculate the correlation function ``C(t)`` for a given bosonic bath and time list.
 
 ## if the input bosonic bath did not apply rotating wave approximation (RWA)
@@ -524,7 +524,7 @@ C^{\nu=\pm}(t)=\sum_i \eta_i^\nu e^{-\gamma_i^\nu t}
 - `cplist::Vector{ComplexF64}` : a list of the value of the absorption (``\nu=+``) correlation function according to the given time list.
 - `cmlist::Vector{ComplexF64}` : a list of the value of the emission (``\nu=-``) correlation function according to the given time list.
 """
-function C(bath::BosonBath, tlist::AbstractVector)
+function correlation_function(bath::BosonBath, tlist::AbstractVector)
     T = (bath[1]).types
 
     # without RWA
@@ -543,7 +543,7 @@ function C(bath::BosonBath, tlist::AbstractVector)
         end
         return clist
 
-        # with RWA
+    # with RWA
     else
         cplist = zeros(ComplexF64, length(tlist))
         cmlist = zeros(ComplexF64, length(tlist))
diff --git a/src/bath/FermionBath.jl b/src/bath/FermionBath.jl
index 7ddd2883..fd2fba6c 100644
--- a/src/bath/FermionBath.jl
+++ b/src/bath/FermionBath.jl
@@ -208,7 +208,7 @@ function fermionEmit(
 end
 
 @doc raw"""
-    C(bath, tlist)
+    correlation_function(bath, tlist)
 Calculate the correlation function ``C^{\nu=+}(t)`` and ``C^{\nu=-}(t)`` for a given fermionic bath and time list.
 Here, ``\nu=+`` represents the absorption process and ``\nu=-`` represents the emission process.
 
@@ -224,7 +224,7 @@ C^{\nu=\pm}(t)=\sum_i \eta_i^\nu e^{-\gamma_i^\nu t}
 - `cplist::Vector{ComplexF64}` : a list of the value of the absorption (``\nu=+``) correlation function according to the given time list.
 - `cmlist::Vector{ComplexF64}` : a list of the value of the emission (``\nu=-``) correlation function according to the given time list.
 """
-function C(bath::FermionBath, tlist::AbstractVector)
+function correlation_function(bath::FermionBath, tlist::AbstractVector)
     cplist = zeros(ComplexF64, length(tlist))
     cmlist = zeros(ComplexF64, length(tlist))
     for (i, t) in enumerate(tlist)
diff --git a/src/deprecated.jl b/src/deprecated.jl
index 2c4ad7db..0d4230d7 100644
--- a/src/deprecated.jl
+++ b/src/deprecated.jl
@@ -14,12 +14,14 @@ function deprecated_foo(args...; kwargs...)
 end
 =#
 
-export evolution, SteadyState
+export evolution, SteadyState, C
 
 SteadyState(args...; kwargs...) = error("`SteadyState` has been deprecated, please use `steadystate` instead.")
 
 evolution(args...; kwargs...) = error("`evolution` has been deprecated, please use `HEOMsolve` instead.")
 
+C(args...; kwargs...) = error("`C` has been deprecated, please use `correlation_function` instead.")
+
 _HEOMSuperOp_deprecated_method_error() = error(
     "Specifying `mul_basis = \"L\", \"R\", or \"LR\"` to `HEOMSuperOp` has been deprecated, try to construct system SuperOperator manually by using `spre`, `spost`, or `sprepost`.",
 )
diff --git a/test/bath.jl b/test/bath.jl
index d12b0fc1..dc77f4cf 100644
--- a/test/bath.jl
+++ b/test/bath.jl
@@ -41,7 +41,8 @@
     b = BosonBath(op, η1, γ1, η2, γ2)
     types = ["bRI", "bRI", "bR", "bR", "bI", "bI", "bI"]
     @test length(b) == 7
-    @test C(b, [0.183183])[1] ≈ C(bs, [0.183183])[1]
+    @test correlation_function(b, 0.183183)[1] ≈ correlation_function(bs, [0.183183])[1]
+    @test_throws ErrorException C(bs, [0.183183])
 
     ## check for η and γ list, and coupling operator
     η = Vector{ComplexF64}(undef, length(b))
@@ -80,7 +81,8 @@
     b = BosonBathRWA(op, η1, γ3, η2, γ4)
     types = ["bA", "bA", "bA", "bA", "bA", "bE", "bE", "bE", "bE", "bE"]
     @test length(bs) == 10
-    cp, cm = C(b, [0.183183])
+    cp, cm = correlation_function(b, 0.183183)
+    @test_throws ErrorException C(b, [0.183183])
     @test cp[1] ≈ 22.384506765987076 + 0.7399082821797519im
     @test cm[1] ≈ 26.994911851482776 - 0.799138487523946im
 
@@ -133,7 +135,8 @@
     b = FermionBath(op, η1, γ3, η2, γ4)
     types = ["fA", "fA", "fA", "fA", "fA", "fE", "fE", "fE", "fE", "fE"]
     @test length(b) == 10
-    cp, cm = C(b, [0.183183])
+    cp, cm = correlation_function(b, 0.183183)
+    @test_throws ErrorException C(b, [0.183183])
     @test cp[1] ≈ 22.384506765987076 + 0.7399082821797519im
     @test cm[1] ≈ 26.994911851482776 - 0.799138487523946im
     for (i, e) in enumerate(b)
@@ -164,7 +167,7 @@
 
     ################################################
     # Exponent
-    @test C(BosonBath(op, η0, γ0), [0.123])[1] ≈ η0[1] * exp(-γ0[1] * 0.123)
+    @test correlation_function(BosonBath(op, η0, γ0), 0.123)[1] ≈ η0[1] * exp(-γ0[1] * 0.123)
 
     ## check exceptions
     @test_throws ErrorException b[11]
diff --git a/test/runtests.jl b/test/runtests.jl
index ef3d1904..5433806f 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -23,12 +23,12 @@ core_tests = [
     "time_evolution.jl",
 ]
 
-if (GROUP == "All") || (GROUP == "Code_Quality")
+#= if (GROUP == "All") || (GROUP == "Code_Quality")
     using HierarchicalEOM
     using Aqua, JET
 
     include(joinpath(testdir, "code_quality.jl"))
-end
+end =#
 
 if (GROUP == "All") || (GROUP == "Core")
     using HierarchicalEOM