diff --git a/datafusion/core/src/dataframe/mod.rs b/datafusion/core/src/dataframe/mod.rs
index 65cf2d700e2c..06a85d303687 100644
--- a/datafusion/core/src/dataframe/mod.rs
+++ b/datafusion/core/src/dataframe/mod.rs
@@ -50,12 +50,12 @@ use datafusion_common::{
 };
 use datafusion_expr::lit;
 use datafusion_expr::{
-    avg, count, max, min, stddev, utils::COUNT_STAR_EXPANSION,
-    TableProviderFilterPushDown, UNNAMED_TABLE,
+    avg, count, max, min, utils::COUNT_STAR_EXPANSION, TableProviderFilterPushDown,
+    UNNAMED_TABLE,
 };
 use datafusion_expr::{case, is_null};
-use datafusion_functions_aggregate::expr_fn::median;
 use datafusion_functions_aggregate::expr_fn::sum;
+use datafusion_functions_aggregate::expr_fn::{median, stddev};
 
 use async_trait::async_trait;
 
diff --git a/datafusion/expr/src/aggregate_function.rs b/datafusion/expr/src/aggregate_function.rs
index 8340fe96ad90..edefd0f3edf2 100644
--- a/datafusion/expr/src/aggregate_function.rs
+++ b/datafusion/expr/src/aggregate_function.rs
@@ -49,10 +49,6 @@ pub enum AggregateFunction {
     NthValue,
     /// Variance (Population)
     VariancePop,
-    /// Standard Deviation (Sample)
-    Stddev,
-    /// Standard Deviation (Population)
-    StddevPop,
     /// Correlation
     Correlation,
     /// Slope from linear regression
@@ -107,8 +103,6 @@ impl AggregateFunction {
             ArrayAgg => "ARRAY_AGG",
             NthValue => "NTH_VALUE",
             VariancePop => "VAR_POP",
-            Stddev => "STDDEV",
-            StddevPop => "STDDEV_POP",
             Correlation => "CORR",
             RegrSlope => "REGR_SLOPE",
             RegrIntercept => "REGR_INTERCEPT",
@@ -159,9 +153,6 @@ impl FromStr for AggregateFunction {
             "string_agg" => AggregateFunction::StringAgg,
             // statistical
             "corr" => AggregateFunction::Correlation,
-            "stddev" => AggregateFunction::Stddev,
-            "stddev_pop" => AggregateFunction::StddevPop,
-            "stddev_samp" => AggregateFunction::Stddev,
             "var_pop" => AggregateFunction::VariancePop,
             "regr_slope" => AggregateFunction::RegrSlope,
             "regr_intercept" => AggregateFunction::RegrIntercept,
@@ -231,8 +222,6 @@ impl AggregateFunction {
             AggregateFunction::Correlation => {
                 correlation_return_type(&coerced_data_types[0])
             }
-            AggregateFunction::Stddev => stddev_return_type(&coerced_data_types[0]),
-            AggregateFunction::StddevPop => stddev_return_type(&coerced_data_types[0]),
             AggregateFunction::RegrSlope
             | AggregateFunction::RegrIntercept
             | AggregateFunction::RegrCount
@@ -304,8 +293,6 @@ impl AggregateFunction {
             }
             AggregateFunction::Avg
             | AggregateFunction::VariancePop
-            | AggregateFunction::Stddev
-            | AggregateFunction::StddevPop
             | AggregateFunction::ApproxMedian => {
                 Signature::uniform(1, NUMERICS.to_vec(), Volatility::Immutable)
             }
diff --git a/datafusion/expr/src/expr_fn.rs b/datafusion/expr/src/expr_fn.rs
index b1d9eb057753..0360478eac54 100644
--- a/datafusion/expr/src/expr_fn.rs
+++ b/datafusion/expr/src/expr_fn.rs
@@ -383,18 +383,6 @@ pub fn scalar_subquery(subquery: Arc<LogicalPlan>) -> Expr {
     })
 }
 
-/// Create an expression to represent the stddev() aggregate function
-pub fn stddev(expr: Expr) -> Expr {
-    Expr::AggregateFunction(AggregateFunction::new(
-        aggregate_function::AggregateFunction::Stddev,
-        vec![expr],
-        false,
-        None,
-        None,
-        None,
-    ))
-}
-
 /// Create a grouping set
 pub fn grouping_set(exprs: Vec<Vec<Expr>>) -> Expr {
     Expr::GroupingSet(GroupingSet::GroupingSets(exprs))
diff --git a/datafusion/expr/src/type_coercion/aggregates.rs b/datafusion/expr/src/type_coercion/aggregates.rs
index 98738bcde7eb..3ac6fa41d9ea 100644
--- a/datafusion/expr/src/type_coercion/aggregates.rs
+++ b/datafusion/expr/src/type_coercion/aggregates.rs
@@ -161,16 +161,6 @@ pub fn coerce_types(
             }
             Ok(vec![Float64, Float64])
         }
-        AggregateFunction::Stddev | AggregateFunction::StddevPop => {
-            if !is_stddev_support_arg_type(&input_types[0]) {
-                return plan_err!(
-                    "The function {:?} does not support inputs of type {:?}.",
-                    agg_fun,
-                    input_types[0]
-                );
-            }
-            Ok(vec![Float64])
-        }
         AggregateFunction::Correlation => {
             if !is_correlation_support_arg_type(&input_types[0]) {
                 return plan_err!(
@@ -408,15 +398,6 @@ pub fn correlation_return_type(arg_type: &DataType) -> Result<DataType> {
     }
 }
 
-/// function return type of standard deviation
-pub fn stddev_return_type(arg_type: &DataType) -> Result<DataType> {
-    if NUMERICS.contains(arg_type) {
-        Ok(DataType::Float64)
-    } else {
-        plan_err!("STDDEV does not support {arg_type:?}")
-    }
-}
-
 /// function return type of an average
 pub fn avg_return_type(arg_type: &DataType) -> Result<DataType> {
     match arg_type {
@@ -511,13 +492,6 @@ pub fn is_covariance_support_arg_type(arg_type: &DataType) -> bool {
     )
 }
 
-pub fn is_stddev_support_arg_type(arg_type: &DataType) -> bool {
-    matches!(
-        arg_type,
-        arg_type if NUMERICS.contains(arg_type)
-    )
-}
-
 pub fn is_correlation_support_arg_type(arg_type: &DataType) -> bool {
     matches!(
         arg_type,
@@ -664,17 +638,6 @@ mod tests {
         Ok(())
     }
 
-    #[test]
-    fn test_stddev_return_data_type() -> Result<()> {
-        let data_type = DataType::Float64;
-        let result_type = stddev_return_type(&data_type)?;
-        assert_eq!(DataType::Float64, result_type);
-
-        let data_type = DataType::Decimal128(36, 10);
-        assert!(stddev_return_type(&data_type).is_err());
-        Ok(())
-    }
-
     #[test]
     fn test_covariance_return_data_type() -> Result<()> {
         let data_type = DataType::Float64;
diff --git a/datafusion/functions-aggregate/src/lib.rs b/datafusion/functions-aggregate/src/lib.rs
index ff02d25ad00b..2f58b9afaccd 100644
--- a/datafusion/functions-aggregate/src/lib.rs
+++ b/datafusion/functions-aggregate/src/lib.rs
@@ -58,6 +58,7 @@ pub mod macros;
 pub mod covariance;
 pub mod first_last;
 pub mod median;
+pub mod stddev;
 pub mod sum;
 pub mod variance;
 
@@ -74,6 +75,8 @@ pub mod expr_fn {
     pub use super::first_last::first_value;
     pub use super::first_last::last_value;
     pub use super::median::median;
+    pub use super::stddev::stddev;
+    pub use super::stddev::stddev_pop;
     pub use super::sum::sum;
     pub use super::variance::var_sample;
 }
@@ -88,6 +91,8 @@ pub fn all_default_aggregate_functions() -> Vec<Arc<AggregateUDF>> {
         covariance::covar_pop_udaf(),
         median::median_udaf(),
         variance::var_samp_udaf(),
+        stddev::stddev_udaf(),
+        stddev::stddev_pop_udaf(),
     ]
 }
 
diff --git a/datafusion/functions-aggregate/src/stddev.rs b/datafusion/functions-aggregate/src/stddev.rs
new file mode 100644
index 000000000000..4c3effe7650a
--- /dev/null
+++ b/datafusion/functions-aggregate/src/stddev.rs
@@ -0,0 +1,366 @@
+// Licensed to the Apache Software Foundation (ASF) under one
+// or more contributor license agreements.  See the NOTICE file
+// distributed with this work for additional information
+// regarding copyright ownership.  The ASF licenses this file
+// to you under the Apache License, Version 2.0 (the
+// "License"); you may not use this file except in compliance
+// with the License.  You may obtain a copy of the License at
+//
+//   http://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing,
+// software distributed under the License is distributed on an
+// "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
+// KIND, either express or implied.  See the License for the
+// specific language governing permissions and limitations
+// under the License.
+
+//! Defines physical expressions that can evaluated at runtime during query execution
+
+use std::any::Any;
+use std::fmt::{Debug, Formatter};
+
+use arrow::{array::ArrayRef, datatypes::DataType, datatypes::Field};
+
+use datafusion_common::{internal_err, not_impl_err, Result};
+use datafusion_common::{plan_err, ScalarValue};
+use datafusion_expr::function::{AccumulatorArgs, StateFieldsArgs};
+use datafusion_expr::utils::format_state_name;
+use datafusion_expr::{Accumulator, AggregateUDFImpl, Signature, Volatility};
+use datafusion_physical_expr_common::aggregate::stats::StatsType;
+
+use crate::variance::VarianceAccumulator;
+
+make_udaf_expr_and_func!(
+    Stddev,
+    stddev,
+    expression,
+    "Compute the standard deviation of a set of numbers",
+    stddev_udaf
+);
+
+/// STDDEV and STDDEV_SAMP (standard deviation) aggregate expression
+pub struct Stddev {
+    signature: Signature,
+    alias: Vec<String>,
+}
+
+impl Debug for Stddev {
+    fn fmt(&self, f: &mut Formatter<'_>) -> std::fmt::Result {
+        f.debug_struct("Stddev")
+            .field("name", &self.name())
+            .field("signature", &self.signature)
+            .finish()
+    }
+}
+
+impl Default for Stddev {
+    fn default() -> Self {
+        Self::new()
+    }
+}
+
+impl Stddev {
+    /// Create a new STDDEV aggregate function
+    pub fn new() -> Self {
+        Self {
+            signature: Signature::numeric(1, Volatility::Immutable),
+            alias: vec!["stddev_samp".to_string()],
+        }
+    }
+}
+
+impl AggregateUDFImpl for Stddev {
+    /// Return a reference to Any that can be used for downcasting
+    fn as_any(&self) -> &dyn Any {
+        self
+    }
+
+    fn name(&self) -> &str {
+        "stddev"
+    }
+
+    fn signature(&self) -> &Signature {
+        &self.signature
+    }
+
+    fn return_type(&self, arg_types: &[DataType]) -> Result<DataType> {
+        if !arg_types[0].is_numeric() {
+            return plan_err!("Stddev requires numeric input types");
+        }
+
+        Ok(DataType::Float64)
+    }
+
+    fn state_fields(&self, args: StateFieldsArgs) -> Result<Vec<Field>> {
+        Ok(vec![
+            Field::new(
+                format_state_name(args.name, "count"),
+                DataType::UInt64,
+                true,
+            ),
+            Field::new(
+                format_state_name(args.name, "mean"),
+                DataType::Float64,
+                true,
+            ),
+            Field::new(format_state_name(args.name, "m2"), DataType::Float64, true),
+        ])
+    }
+
+    fn accumulator(&self, acc_args: AccumulatorArgs) -> Result<Box<dyn Accumulator>> {
+        if acc_args.is_distinct {
+            return not_impl_err!("STDDEV_POP(DISTINCT) aggregations are not available");
+        }
+        Ok(Box::new(StddevAccumulator::try_new(StatsType::Sample)?))
+    }
+
+    fn aliases(&self) -> &[String] {
+        &self.alias
+    }
+}
+
+make_udaf_expr_and_func!(
+    StddevPop,
+    stddev_pop,
+    expression,
+    "Compute the population standard deviation of a set of numbers",
+    stddev_pop_udaf
+);
+
+/// STDDEV_POP population aggregate expression
+pub struct StddevPop {
+    signature: Signature,
+}
+
+impl Debug for StddevPop {
+    fn fmt(&self, f: &mut Formatter<'_>) -> std::fmt::Result {
+        f.debug_struct("StddevPop")
+            .field("name", &self.name())
+            .field("signature", &self.signature)
+            .finish()
+    }
+}
+
+impl Default for StddevPop {
+    fn default() -> Self {
+        Self::new()
+    }
+}
+
+impl StddevPop {
+    /// Create a new STDDEV_POP aggregate function
+    pub fn new() -> Self {
+        Self {
+            signature: Signature::numeric(1, Volatility::Immutable),
+        }
+    }
+}
+
+impl AggregateUDFImpl for StddevPop {
+    /// Return a reference to Any that can be used for downcasting
+    fn as_any(&self) -> &dyn Any {
+        self
+    }
+
+    fn name(&self) -> &str {
+        "stddev_pop"
+    }
+
+    fn signature(&self) -> &Signature {
+        &self.signature
+    }
+
+    fn state_fields(&self, args: StateFieldsArgs) -> Result<Vec<Field>> {
+        Ok(vec![
+            Field::new(
+                format_state_name(args.name, "count"),
+                DataType::UInt64,
+                true,
+            ),
+            Field::new(
+                format_state_name(args.name, "mean"),
+                DataType::Float64,
+                true,
+            ),
+            Field::new(format_state_name(args.name, "m2"), DataType::Float64, true),
+        ])
+    }
+
+    fn accumulator(&self, acc_args: AccumulatorArgs) -> Result<Box<dyn Accumulator>> {
+        if acc_args.is_distinct {
+            return not_impl_err!("STDDEV_POP(DISTINCT) aggregations are not available");
+        }
+        Ok(Box::new(StddevAccumulator::try_new(StatsType::Population)?))
+    }
+
+    fn return_type(&self, arg_types: &[DataType]) -> Result<DataType> {
+        if !arg_types[0].is_numeric() {
+            return plan_err!("StddevPop requires numeric input types");
+        }
+
+        Ok(DataType::Float64)
+    }
+}
+
+/// An accumulator to compute the average
+#[derive(Debug)]
+pub struct StddevAccumulator {
+    variance: VarianceAccumulator,
+}
+
+impl StddevAccumulator {
+    /// Creates a new `StddevAccumulator`
+    pub fn try_new(s_type: StatsType) -> Result<Self> {
+        Ok(Self {
+            variance: VarianceAccumulator::try_new(s_type)?,
+        })
+    }
+
+    pub fn get_m2(&self) -> f64 {
+        self.variance.get_m2()
+    }
+}
+
+impl Accumulator for StddevAccumulator {
+    fn state(&mut self) -> Result<Vec<ScalarValue>> {
+        Ok(vec![
+            ScalarValue::from(self.variance.get_count()),
+            ScalarValue::from(self.variance.get_mean()),
+            ScalarValue::from(self.variance.get_m2()),
+        ])
+    }
+
+    fn update_batch(&mut self, values: &[ArrayRef]) -> Result<()> {
+        self.variance.update_batch(values)
+    }
+
+    fn retract_batch(&mut self, values: &[ArrayRef]) -> Result<()> {
+        self.variance.retract_batch(values)
+    }
+
+    fn merge_batch(&mut self, states: &[ArrayRef]) -> Result<()> {
+        self.variance.merge_batch(states)
+    }
+
+    fn evaluate(&mut self) -> Result<ScalarValue> {
+        let variance = self.variance.evaluate()?;
+        match variance {
+            ScalarValue::Float64(e) => {
+                if e.is_none() {
+                    Ok(ScalarValue::Float64(None))
+                } else {
+                    Ok(ScalarValue::Float64(e.map(|f| f.sqrt())))
+                }
+            }
+            _ => internal_err!("Variance should be f64"),
+        }
+    }
+
+    fn size(&self) -> usize {
+        std::mem::align_of_val(self) - std::mem::align_of_val(&self.variance)
+            + self.variance.size()
+    }
+
+    fn supports_retract_batch(&self) -> bool {
+        self.variance.supports_retract_batch()
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use std::sync::Arc;
+
+    use arrow::{array::*, datatypes::*};
+
+    use datafusion_expr::AggregateUDF;
+    use datafusion_physical_expr_common::aggregate::utils::get_accum_scalar_values_as_arrays;
+    use datafusion_physical_expr_common::expressions::column::col;
+
+    use super::*;
+
+    #[test]
+    fn stddev_f64_merge_1() -> Result<()> {
+        let a = Arc::new(Float64Array::from(vec![1_f64, 2_f64, 3_f64]));
+        let b = Arc::new(Float64Array::from(vec![4_f64, 5_f64]));
+
+        let schema = Schema::new(vec![Field::new("a", DataType::Float64, false)]);
+
+        let batch1 = RecordBatch::try_new(Arc::new(schema.clone()), vec![a])?;
+        let batch2 = RecordBatch::try_new(Arc::new(schema.clone()), vec![b])?;
+
+        let agg1 = stddev_pop_udaf();
+        let agg2 = stddev_pop_udaf();
+
+        let actual = merge(&batch1, &batch2, agg1, agg2, &schema)?;
+        assert_eq!(actual, ScalarValue::from(std::f64::consts::SQRT_2));
+
+        Ok(())
+    }
+
+    #[test]
+    fn stddev_f64_merge_2() -> Result<()> {
+        let a = Arc::new(Float64Array::from(vec![1_f64, 2_f64, 3_f64, 4_f64, 5_f64]));
+        let b = Arc::new(Float64Array::from(vec![None]));
+
+        let schema = Schema::new(vec![Field::new("a", DataType::Float64, true)]);
+
+        let batch1 = RecordBatch::try_new(Arc::new(schema.clone()), vec![a])?;
+        let batch2 = RecordBatch::try_new(Arc::new(schema.clone()), vec![b])?;
+
+        let agg1 = stddev_pop_udaf();
+        let agg2 = stddev_pop_udaf();
+
+        let actual = merge(&batch1, &batch2, agg1, agg2, &schema)?;
+        assert_eq!(actual, ScalarValue::from(std::f64::consts::SQRT_2));
+
+        Ok(())
+    }
+
+    fn merge(
+        batch1: &RecordBatch,
+        batch2: &RecordBatch,
+        agg1: Arc<AggregateUDF>,
+        agg2: Arc<AggregateUDF>,
+        schema: &Schema,
+    ) -> Result<ScalarValue> {
+        let args1 = AccumulatorArgs {
+            data_type: &DataType::Float64,
+            schema,
+            ignore_nulls: false,
+            sort_exprs: &[],
+            name: "a",
+            is_distinct: false,
+            input_type: &DataType::Float64,
+            args_num: 1,
+        };
+
+        let args2 = AccumulatorArgs {
+            data_type: &DataType::Float64,
+            schema,
+            ignore_nulls: false,
+            sort_exprs: &[],
+            name: "a",
+            is_distinct: false,
+            input_type: &DataType::Float64,
+            args_num: 1,
+        };
+
+        let mut accum1 = agg1.accumulator(args1)?;
+        let mut accum2 = agg2.accumulator(args2)?;
+
+        let value1 = vec![col("a", schema)?
+            .evaluate(batch1)
+            .and_then(|v| v.into_array(batch1.num_rows()))?];
+        let value2 = vec![col("a", schema)?
+            .evaluate(batch2)
+            .and_then(|v| v.into_array(batch2.num_rows()))?];
+
+        accum1.update_batch(&value1)?;
+        accum2.update_batch(&value2)?;
+        let state2 = get_accum_scalar_values_as_arrays(accum2.as_mut())?;
+        accum1.merge_batch(&state2)?;
+        let result = accum1.evaluate()?;
+        Ok(result)
+    }
+}
diff --git a/datafusion/physical-expr/src/aggregate/build_in.rs b/datafusion/physical-expr/src/aggregate/build_in.rs
index 9c40dc4520c4..d92d0cd61c65 100644
--- a/datafusion/physical-expr/src/aggregate/build_in.rs
+++ b/datafusion/physical-expr/src/aggregate/build_in.rs
@@ -163,22 +163,6 @@ pub fn create_aggregate_expr(
         (AggregateFunction::VariancePop, true) => {
             return not_impl_err!("VAR_POP(DISTINCT) aggregations are not available");
         }
-        (AggregateFunction::Stddev, false) => Arc::new(expressions::Stddev::new(
-            input_phy_exprs[0].clone(),
-            name,
-            data_type,
-        )),
-        (AggregateFunction::Stddev, true) => {
-            return not_impl_err!("STDDEV(DISTINCT) aggregations are not available");
-        }
-        (AggregateFunction::StddevPop, false) => Arc::new(expressions::StddevPop::new(
-            input_phy_exprs[0].clone(),
-            name,
-            data_type,
-        )),
-        (AggregateFunction::StddevPop, true) => {
-            return not_impl_err!("STDDEV_POP(DISTINCT) aggregations are not available");
-        }
         (AggregateFunction::Correlation, false) => {
             Arc::new(expressions::Correlation::new(
                 input_phy_exprs[0].clone(),
@@ -356,13 +340,13 @@ pub fn create_aggregate_expr(
 #[cfg(test)]
 mod tests {
     use arrow::datatypes::{DataType, Field};
-    use expressions::{StddevPop, VariancePop};
+    use expressions::VariancePop;
 
     use super::*;
     use crate::expressions::{
         try_cast, ApproxDistinct, ApproxMedian, ApproxPercentileCont, ArrayAgg, Avg,
         BitAnd, BitOr, BitXor, BoolAnd, BoolOr, Count, DistinctArrayAgg, DistinctCount,
-        Max, Min, Stddev,
+        Max, Min,
     };
 
     use datafusion_common::{plan_err, DataFusionError, ScalarValue};
@@ -747,82 +731,6 @@ mod tests {
         Ok(())
     }
 
-    #[test]
-    fn test_stddev_expr() -> Result<()> {
-        let funcs = vec![AggregateFunction::Stddev];
-        let data_types = vec![
-            DataType::UInt32,
-            DataType::UInt64,
-            DataType::Int32,
-            DataType::Int64,
-            DataType::Float32,
-            DataType::Float64,
-        ];
-        for fun in funcs {
-            for data_type in &data_types {
-                let input_schema =
-                    Schema::new(vec![Field::new("c1", data_type.clone(), true)]);
-                let input_phy_exprs: Vec<Arc<dyn PhysicalExpr>> = vec![Arc::new(
-                    expressions::Column::new_with_schema("c1", &input_schema).unwrap(),
-                )];
-                let result_agg_phy_exprs = create_physical_agg_expr_for_test(
-                    &fun,
-                    false,
-                    &input_phy_exprs[0..1],
-                    &input_schema,
-                    "c1",
-                )?;
-                if fun == AggregateFunction::Stddev {
-                    assert!(result_agg_phy_exprs.as_any().is::<Stddev>());
-                    assert_eq!("c1", result_agg_phy_exprs.name());
-                    assert_eq!(
-                        Field::new("c1", DataType::Float64, true),
-                        result_agg_phy_exprs.field().unwrap()
-                    )
-                }
-            }
-        }
-        Ok(())
-    }
-
-    #[test]
-    fn test_stddev_pop_expr() -> Result<()> {
-        let funcs = vec![AggregateFunction::StddevPop];
-        let data_types = vec![
-            DataType::UInt32,
-            DataType::UInt64,
-            DataType::Int32,
-            DataType::Int64,
-            DataType::Float32,
-            DataType::Float64,
-        ];
-        for fun in funcs {
-            for data_type in &data_types {
-                let input_schema =
-                    Schema::new(vec![Field::new("c1", data_type.clone(), true)]);
-                let input_phy_exprs: Vec<Arc<dyn PhysicalExpr>> = vec![Arc::new(
-                    expressions::Column::new_with_schema("c1", &input_schema).unwrap(),
-                )];
-                let result_agg_phy_exprs = create_physical_agg_expr_for_test(
-                    &fun,
-                    false,
-                    &input_phy_exprs[0..1],
-                    &input_schema,
-                    "c1",
-                )?;
-                if fun == AggregateFunction::StddevPop {
-                    assert!(result_agg_phy_exprs.as_any().is::<StddevPop>());
-                    assert_eq!("c1", result_agg_phy_exprs.name());
-                    assert_eq!(
-                        Field::new("c1", DataType::Float64, true),
-                        result_agg_phy_exprs.field().unwrap()
-                    )
-                }
-            }
-        }
-        Ok(())
-    }
-
     #[test]
     fn test_median_expr() -> Result<()> {
         let funcs = vec![AggregateFunction::ApproxMedian];
@@ -939,32 +847,6 @@ mod tests {
         assert!(observed.is_err());
     }
 
-    #[test]
-    fn test_stddev_return_type() -> Result<()> {
-        let observed = AggregateFunction::Stddev.return_type(&[DataType::Float32])?;
-        assert_eq!(DataType::Float64, observed);
-
-        let observed = AggregateFunction::Stddev.return_type(&[DataType::Float64])?;
-        assert_eq!(DataType::Float64, observed);
-
-        let observed = AggregateFunction::Stddev.return_type(&[DataType::Int32])?;
-        assert_eq!(DataType::Float64, observed);
-
-        let observed = AggregateFunction::Stddev.return_type(&[DataType::UInt32])?;
-        assert_eq!(DataType::Float64, observed);
-
-        let observed = AggregateFunction::Stddev.return_type(&[DataType::Int64])?;
-        assert_eq!(DataType::Float64, observed);
-
-        Ok(())
-    }
-
-    #[test]
-    fn test_stddev_no_utf8() {
-        let observed = AggregateFunction::Stddev.return_type(&[DataType::Utf8]);
-        assert!(observed.is_err());
-    }
-
     // Helper function
     // Create aggregate expr with type coercion
     fn create_physical_agg_expr_for_test(
diff --git a/datafusion/physical-expr/src/aggregate/stddev.rs b/datafusion/physical-expr/src/aggregate/stddev.rs
index ec8d8cea67c4..3ade67b51905 100644
--- a/datafusion/physical-expr/src/aggregate/stddev.rs
+++ b/datafusion/physical-expr/src/aggregate/stddev.rs
@@ -17,168 +17,14 @@
 
 //! Defines physical expressions that can evaluated at runtime during query execution
 
-use std::any::Any;
-use std::sync::Arc;
+use arrow::array::ArrayRef;
 
-use crate::aggregate::stats::StatsType;
-use crate::aggregate::utils::down_cast_any_ref;
-use crate::aggregate::variance::VarianceAccumulator;
-use crate::expressions::format_state_name;
-use crate::{AggregateExpr, PhysicalExpr};
-use arrow::{array::ArrayRef, datatypes::DataType, datatypes::Field};
 use datafusion_common::ScalarValue;
 use datafusion_common::{internal_err, Result};
 use datafusion_expr::Accumulator;
 
-/// STDDEV and STDDEV_SAMP (standard deviation) aggregate expression
-#[derive(Debug)]
-pub struct Stddev {
-    name: String,
-    expr: Arc<dyn PhysicalExpr>,
-}
-
-/// STDDEV_POP population aggregate expression
-#[derive(Debug)]
-pub struct StddevPop {
-    name: String,
-    expr: Arc<dyn PhysicalExpr>,
-}
-
-impl Stddev {
-    /// Create a new STDDEV aggregate function
-    pub fn new(
-        expr: Arc<dyn PhysicalExpr>,
-        name: impl Into<String>,
-        data_type: DataType,
-    ) -> Self {
-        // the result of stddev just support FLOAT64 and Decimal data type.
-        assert!(matches!(data_type, DataType::Float64));
-        Self {
-            name: name.into(),
-            expr,
-        }
-    }
-}
-
-impl AggregateExpr for Stddev {
-    /// Return a reference to Any that can be used for downcasting
-    fn as_any(&self) -> &dyn Any {
-        self
-    }
-
-    fn field(&self) -> Result<Field> {
-        Ok(Field::new(&self.name, DataType::Float64, true))
-    }
-
-    fn create_accumulator(&self) -> Result<Box<dyn Accumulator>> {
-        Ok(Box::new(StddevAccumulator::try_new(StatsType::Sample)?))
-    }
-
-    fn create_sliding_accumulator(&self) -> Result<Box<dyn Accumulator>> {
-        Ok(Box::new(StddevAccumulator::try_new(StatsType::Sample)?))
-    }
-
-    fn state_fields(&self) -> Result<Vec<Field>> {
-        Ok(vec![
-            Field::new(
-                format_state_name(&self.name, "count"),
-                DataType::UInt64,
-                true,
-            ),
-            Field::new(
-                format_state_name(&self.name, "mean"),
-                DataType::Float64,
-                true,
-            ),
-            Field::new(format_state_name(&self.name, "m2"), DataType::Float64, true),
-        ])
-    }
-
-    fn expressions(&self) -> Vec<Arc<dyn PhysicalExpr>> {
-        vec![self.expr.clone()]
-    }
-
-    fn name(&self) -> &str {
-        &self.name
-    }
-}
-
-impl PartialEq<dyn Any> for Stddev {
-    fn eq(&self, other: &dyn Any) -> bool {
-        down_cast_any_ref(other)
-            .downcast_ref::<Self>()
-            .map(|x| self.name == x.name && self.expr.eq(&x.expr))
-            .unwrap_or(false)
-    }
-}
-
-impl StddevPop {
-    /// Create a new STDDEV aggregate function
-    pub fn new(
-        expr: Arc<dyn PhysicalExpr>,
-        name: impl Into<String>,
-        data_type: DataType,
-    ) -> Self {
-        // the result of stddev just support FLOAT64 and Decimal data type.
-        assert!(matches!(data_type, DataType::Float64));
-        Self {
-            name: name.into(),
-            expr,
-        }
-    }
-}
-
-impl AggregateExpr for StddevPop {
-    /// Return a reference to Any that can be used for downcasting
-    fn as_any(&self) -> &dyn Any {
-        self
-    }
-
-    fn field(&self) -> Result<Field> {
-        Ok(Field::new(&self.name, DataType::Float64, true))
-    }
-
-    fn create_accumulator(&self) -> Result<Box<dyn Accumulator>> {
-        Ok(Box::new(StddevAccumulator::try_new(StatsType::Population)?))
-    }
-
-    fn create_sliding_accumulator(&self) -> Result<Box<dyn Accumulator>> {
-        Ok(Box::new(StddevAccumulator::try_new(StatsType::Population)?))
-    }
-
-    fn state_fields(&self) -> Result<Vec<Field>> {
-        Ok(vec![
-            Field::new(
-                format_state_name(&self.name, "count"),
-                DataType::UInt64,
-                true,
-            ),
-            Field::new(
-                format_state_name(&self.name, "mean"),
-                DataType::Float64,
-                true,
-            ),
-            Field::new(format_state_name(&self.name, "m2"), DataType::Float64, true),
-        ])
-    }
-
-    fn expressions(&self) -> Vec<Arc<dyn PhysicalExpr>> {
-        vec![self.expr.clone()]
-    }
-
-    fn name(&self) -> &str {
-        &self.name
-    }
-}
-
-impl PartialEq<dyn Any> for StddevPop {
-    fn eq(&self, other: &dyn Any) -> bool {
-        down_cast_any_ref(other)
-            .downcast_ref::<Self>()
-            .map(|x| self.name == x.name && self.expr.eq(&x.expr))
-            .unwrap_or(false)
-    }
-}
+use crate::aggregate::stats::StatsType;
+use crate::aggregate::variance::VarianceAccumulator;
 
 /// An accumulator to compute the average
 #[derive(Debug)]
@@ -239,99 +85,3 @@ impl Accumulator for StddevAccumulator {
             + self.variance.size()
     }
 }
-
-#[cfg(test)]
-mod tests {
-    use super::*;
-    use crate::aggregate::utils::get_accum_scalar_values_as_arrays;
-    use crate::expressions::col;
-    use arrow::{array::*, datatypes::*};
-
-    #[test]
-    fn stddev_f64_merge_1() -> Result<()> {
-        let a = Arc::new(Float64Array::from(vec![1_f64, 2_f64, 3_f64]));
-        let b = Arc::new(Float64Array::from(vec![4_f64, 5_f64]));
-
-        let schema = Schema::new(vec![Field::new("a", DataType::Float64, false)]);
-
-        let batch1 = RecordBatch::try_new(Arc::new(schema.clone()), vec![a])?;
-        let batch2 = RecordBatch::try_new(Arc::new(schema.clone()), vec![b])?;
-
-        let agg1 = Arc::new(StddevPop::new(
-            col("a", &schema)?,
-            "bla".to_string(),
-            DataType::Float64,
-        ));
-
-        let agg2 = Arc::new(StddevPop::new(
-            col("a", &schema)?,
-            "bla".to_string(),
-            DataType::Float64,
-        ));
-
-        let actual = merge(&batch1, &batch2, agg1, agg2)?;
-        assert!(actual == ScalarValue::from(std::f64::consts::SQRT_2));
-
-        Ok(())
-    }
-
-    #[test]
-    fn stddev_f64_merge_2() -> Result<()> {
-        let a = Arc::new(Float64Array::from(vec![1_f64, 2_f64, 3_f64, 4_f64, 5_f64]));
-        let b = Arc::new(Float64Array::from(vec![None]));
-
-        let schema = Schema::new(vec![Field::new("a", DataType::Float64, true)]);
-
-        let batch1 = RecordBatch::try_new(Arc::new(schema.clone()), vec![a])?;
-        let batch2 = RecordBatch::try_new(Arc::new(schema.clone()), vec![b])?;
-
-        let agg1 = Arc::new(StddevPop::new(
-            col("a", &schema)?,
-            "bla".to_string(),
-            DataType::Float64,
-        ));
-
-        let agg2 = Arc::new(StddevPop::new(
-            col("a", &schema)?,
-            "bla".to_string(),
-            DataType::Float64,
-        ));
-
-        let actual = merge(&batch1, &batch2, agg1, agg2)?;
-        assert!(actual == ScalarValue::from(std::f64::consts::SQRT_2));
-
-        Ok(())
-    }
-
-    fn merge(
-        batch1: &RecordBatch,
-        batch2: &RecordBatch,
-        agg1: Arc<dyn AggregateExpr>,
-        agg2: Arc<dyn AggregateExpr>,
-    ) -> Result<ScalarValue> {
-        let mut accum1 = agg1.create_accumulator()?;
-        let mut accum2 = agg2.create_accumulator()?;
-        let expr1 = agg1.expressions();
-        let expr2 = agg2.expressions();
-
-        let values1 = expr1
-            .iter()
-            .map(|e| {
-                e.evaluate(batch1)
-                    .and_then(|v| v.into_array(batch1.num_rows()))
-            })
-            .collect::<Result<Vec<_>>>()?;
-        let values2 = expr2
-            .iter()
-            .map(|e| {
-                e.evaluate(batch2)
-                    .and_then(|v| v.into_array(batch2.num_rows()))
-            })
-            .collect::<Result<Vec<_>>>()?;
-        accum1.update_batch(&values1)?;
-        accum2.update_batch(&values2)?;
-        let state2 = get_accum_scalar_values_as_arrays(accum2.as_mut())?;
-        accum1.merge_batch(&state2)?;
-        accum1.evaluate()
-    }
-}
diff --git a/datafusion/physical-expr/src/expressions/mod.rs b/datafusion/physical-expr/src/expressions/mod.rs
index eb53c49d06ce..a6133eeb2511 100644
--- a/datafusion/physical-expr/src/expressions/mod.rs
+++ b/datafusion/physical-expr/src/expressions/mod.rs
@@ -56,7 +56,6 @@ pub use crate::aggregate::min_max::{Max, MaxAccumulator, Min, MinAccumulator};
 pub use crate::aggregate::nth_value::NthValueAgg;
 pub use crate::aggregate::regr::{Regr, RegrType};
 pub use crate::aggregate::stats::StatsType;
-pub use crate::aggregate::stddev::{Stddev, StddevPop};
 pub use crate::aggregate::string_agg::StringAgg;
 pub use crate::aggregate::variance::VariancePop;
 pub use crate::window::cume_dist::{cume_dist, CumeDist};
diff --git a/datafusion/proto/proto/datafusion.proto b/datafusion/proto/proto/datafusion.proto
index 2cda638bc835..5aaae07f4d68 100644
--- a/datafusion/proto/proto/datafusion.proto
+++ b/datafusion/proto/proto/datafusion.proto
@@ -483,8 +483,8 @@ enum AggregateFunction {
   VARIANCE_POP = 8;
   // COVARIANCE = 9;
   // COVARIANCE_POP = 10;
-  STDDEV = 11;
-  STDDEV_POP = 12;
+  // STDDEV = 11;
+  // STDDEV_POP = 12;
   CORRELATION = 13;
   APPROX_PERCENTILE_CONT = 14;
   APPROX_MEDIAN = 15;
diff --git a/datafusion/proto/src/generated/pbjson.rs b/datafusion/proto/src/generated/pbjson.rs
index 077ba5b920cc..cd754e4d9f27 100644
--- a/datafusion/proto/src/generated/pbjson.rs
+++ b/datafusion/proto/src/generated/pbjson.rs
@@ -539,8 +539,6 @@ impl serde::Serialize for AggregateFunction {
             Self::ApproxDistinct => "APPROX_DISTINCT",
             Self::ArrayAgg => "ARRAY_AGG",
             Self::VariancePop => "VARIANCE_POP",
-            Self::Stddev => "STDDEV",
-            Self::StddevPop => "STDDEV_POP",
             Self::Correlation => "CORRELATION",
             Self::ApproxPercentileCont => "APPROX_PERCENTILE_CONT",
             Self::ApproxMedian => "APPROX_MEDIAN",
@@ -580,8 +578,6 @@ impl<'de> serde::Deserialize<'de> for AggregateFunction {
             "APPROX_DISTINCT",
             "ARRAY_AGG",
             "VARIANCE_POP",
-            "STDDEV",
-            "STDDEV_POP",
             "CORRELATION",
             "APPROX_PERCENTILE_CONT",
             "APPROX_MEDIAN",
@@ -650,8 +646,6 @@ impl<'de> serde::Deserialize<'de> for AggregateFunction {
                     "APPROX_DISTINCT" => Ok(AggregateFunction::ApproxDistinct),
                     "ARRAY_AGG" => Ok(AggregateFunction::ArrayAgg),
                     "VARIANCE_POP" => Ok(AggregateFunction::VariancePop),
-                    "STDDEV" => Ok(AggregateFunction::Stddev),
-                    "STDDEV_POP" => Ok(AggregateFunction::StddevPop),
                     "CORRELATION" => Ok(AggregateFunction::Correlation),
                     "APPROX_PERCENTILE_CONT" => Ok(AggregateFunction::ApproxPercentileCont),
                     "APPROX_MEDIAN" => Ok(AggregateFunction::ApproxMedian),
diff --git a/datafusion/proto/src/generated/prost.rs b/datafusion/proto/src/generated/prost.rs
index 76e3e17c7c30..1b38168ba1d2 100644
--- a/datafusion/proto/src/generated/prost.rs
+++ b/datafusion/proto/src/generated/prost.rs
@@ -1927,8 +1927,8 @@ pub enum AggregateFunction {
     VariancePop = 8,
     /// COVARIANCE = 9;
     /// COVARIANCE_POP = 10;
-    Stddev = 11,
-    StddevPop = 12,
+    /// STDDEV = 11;
+    /// STDDEV_POP = 12;
     Correlation = 13,
     ApproxPercentileCont = 14,
     ApproxMedian = 15,
@@ -1966,8 +1966,6 @@ impl AggregateFunction {
             AggregateFunction::ApproxDistinct => "APPROX_DISTINCT",
             AggregateFunction::ArrayAgg => "ARRAY_AGG",
             AggregateFunction::VariancePop => "VARIANCE_POP",
-            AggregateFunction::Stddev => "STDDEV",
-            AggregateFunction::StddevPop => "STDDEV_POP",
             AggregateFunction::Correlation => "CORRELATION",
             AggregateFunction::ApproxPercentileCont => "APPROX_PERCENTILE_CONT",
             AggregateFunction::ApproxMedian => "APPROX_MEDIAN",
@@ -2003,8 +2001,6 @@ impl AggregateFunction {
             "APPROX_DISTINCT" => Some(Self::ApproxDistinct),
             "ARRAY_AGG" => Some(Self::ArrayAgg),
             "VARIANCE_POP" => Some(Self::VariancePop),
-            "STDDEV" => Some(Self::Stddev),
-            "STDDEV_POP" => Some(Self::StddevPop),
             "CORRELATION" => Some(Self::Correlation),
             "APPROX_PERCENTILE_CONT" => Some(Self::ApproxPercentileCont),
             "APPROX_MEDIAN" => Some(Self::ApproxMedian),
diff --git a/datafusion/proto/src/logical_plan/from_proto.rs b/datafusion/proto/src/logical_plan/from_proto.rs
index 0202124f2ef4..0ebf3f71174f 100644
--- a/datafusion/proto/src/logical_plan/from_proto.rs
+++ b/datafusion/proto/src/logical_plan/from_proto.rs
@@ -149,8 +149,6 @@ impl From<protobuf::AggregateFunction> for AggregateFunction {
             protobuf::AggregateFunction::ApproxDistinct => Self::ApproxDistinct,
             protobuf::AggregateFunction::ArrayAgg => Self::ArrayAgg,
             protobuf::AggregateFunction::VariancePop => Self::VariancePop,
-            protobuf::AggregateFunction::Stddev => Self::Stddev,
-            protobuf::AggregateFunction::StddevPop => Self::StddevPop,
             protobuf::AggregateFunction::Correlation => Self::Correlation,
             protobuf::AggregateFunction::RegrSlope => Self::RegrSlope,
             protobuf::AggregateFunction::RegrIntercept => Self::RegrIntercept,
diff --git a/datafusion/proto/src/logical_plan/to_proto.rs b/datafusion/proto/src/logical_plan/to_proto.rs
index 279846c77e32..490ef9a4f240 100644
--- a/datafusion/proto/src/logical_plan/to_proto.rs
+++ b/datafusion/proto/src/logical_plan/to_proto.rs
@@ -120,8 +120,6 @@ impl From<&AggregateFunction> for protobuf::AggregateFunction {
             AggregateFunction::ApproxDistinct => Self::ApproxDistinct,
             AggregateFunction::ArrayAgg => Self::ArrayAgg,
             AggregateFunction::VariancePop => Self::VariancePop,
-            AggregateFunction::Stddev => Self::Stddev,
-            AggregateFunction::StddevPop => Self::StddevPop,
             AggregateFunction::Correlation => Self::Correlation,
             AggregateFunction::RegrSlope => Self::RegrSlope,
             AggregateFunction::RegrIntercept => Self::RegrIntercept,
@@ -418,10 +416,6 @@ pub fn serialize_expr(
                     AggregateFunction::VariancePop => {
                         protobuf::AggregateFunction::VariancePop
                     }
-                    AggregateFunction::Stddev => protobuf::AggregateFunction::Stddev,
-                    AggregateFunction::StddevPop => {
-                        protobuf::AggregateFunction::StddevPop
-                    }
                     AggregateFunction::Correlation => {
                         protobuf::AggregateFunction::Correlation
                     }
diff --git a/datafusion/proto/src/physical_plan/to_proto.rs b/datafusion/proto/src/physical_plan/to_proto.rs
index dd361904592c..66405d4b9af6 100644
--- a/datafusion/proto/src/physical_plan/to_proto.rs
+++ b/datafusion/proto/src/physical_plan/to_proto.rs
@@ -28,8 +28,8 @@ use datafusion::physical_plan::expressions::{
     CastExpr, Column, Correlation, Count, CumeDist, DistinctArrayAgg, DistinctBitXor,
     DistinctCount, Grouping, InListExpr, IsNotNullExpr, IsNullExpr, Literal, Max, Min,
     NegativeExpr, NotExpr, NthValue, NthValueAgg, Ntile, OrderSensitiveArrayAgg, Rank,
-    RankType, Regr, RegrType, RowNumber, Stddev, StddevPop, StringAgg, TryCastExpr,
-    VariancePop, WindowShift,
+    RankType, Regr, RegrType, RowNumber, StringAgg, TryCastExpr, VariancePop,
+    WindowShift,
 };
 use datafusion::physical_plan::udaf::AggregateFunctionExpr;
 use datafusion::physical_plan::windows::{BuiltInWindowExpr, PlainAggregateWindowExpr};
@@ -278,10 +278,6 @@ fn aggr_expr_to_aggr_fn(expr: &dyn AggregateExpr) -> Result<AggrFn> {
         protobuf::AggregateFunction::Avg
     } else if aggr_expr.downcast_ref::<VariancePop>().is_some() {
         protobuf::AggregateFunction::VariancePop
-    } else if aggr_expr.downcast_ref::<Stddev>().is_some() {
-        protobuf::AggregateFunction::Stddev
-    } else if aggr_expr.downcast_ref::<StddevPop>().is_some() {
-        protobuf::AggregateFunction::StddevPop
     } else if aggr_expr.downcast_ref::<Correlation>().is_some() {
         protobuf::AggregateFunction::Correlation
     } else if let Some(regr_expr) = aggr_expr.downcast_ref::<Regr>() {
diff --git a/datafusion/proto/tests/cases/roundtrip_logical_plan.rs b/datafusion/proto/tests/cases/roundtrip_logical_plan.rs
index f32f4b04938f..d0e08033727d 100644
--- a/datafusion/proto/tests/cases/roundtrip_logical_plan.rs
+++ b/datafusion/proto/tests/cases/roundtrip_logical_plan.rs
@@ -32,7 +32,7 @@ use datafusion::execution::context::SessionState;
 use datafusion::execution::runtime_env::{RuntimeConfig, RuntimeEnv};
 use datafusion::execution::FunctionRegistry;
 use datafusion::functions_aggregate::expr_fn::{
-    covar_pop, covar_samp, first_value, median, sum, var_sample,
+    covar_pop, covar_samp, first_value, median, stddev, stddev_pop, sum, var_sample,
 };
 use datafusion::prelude::*;
 use datafusion::test_util::{TestTableFactory, TestTableProvider};
@@ -653,6 +653,8 @@ async fn roundtrip_expr_api() -> Result<()> {
         sum(lit(1)),
         median(lit(2)),
         var_sample(lit(2.2)),
+        stddev(lit(2.2)),
+        stddev_pop(lit(2.2)),
     ];
 
     // ensure expressions created with the expr api can be round tripped
diff --git a/datafusion/sqllogictest/test_files/functions.slt b/datafusion/sqllogictest/test_files/functions.slt
index 56f30ad36296..f04d76822124 100644
--- a/datafusion/sqllogictest/test_files/functions.slt
+++ b/datafusion/sqllogictest/test_files/functions.slt
@@ -491,7 +491,7 @@ statement error Did you mean 'COUNT'?
 SELECT counter(*) from test;
 
 # Aggregate function
-statement error Did you mean 'STDDEV'?
+statement error Did you mean 'stddev'?
 SELECT STDEV(v1) from test;
 
 # Aggregate function
