diff --git a/pandas/core/col.py b/pandas/core/col.py
index 39c4a7fd016c2..b8e965845b4b6 100644
--- a/pandas/core/col.py
+++ b/pandas/core/col.py
@@ -37,6 +37,12 @@
     "__lt__": "<",
     "__eq__": "==",
     "__ne__": "!=",
+    "__and__": "&",
+    "__rand__": "&",
+    "__or__": "|",
+    "__ror__": "|",
+    "__xor__": "^",
+    "__rxor__": "^",
 }
 
 
@@ -157,6 +163,28 @@ def __mod__(self, other: Any) -> Expression:
     def __rmod__(self, other: Any) -> Expression:
         return self._with_binary_op("__rmod__", other)
 
+    # Logical ops
+    def __and__(self, other: Any) -> Expression:
+        return self._with_binary_op("__and__", other)
+
+    def __rand__(self, other: Any) -> Expression:
+        return self._with_binary_op("__rand__", other)
+
+    def __or__(self, other: Any) -> Expression:
+        return self._with_binary_op("__or__", other)
+
+    def __ror__(self, other: Any) -> Expression:
+        return self._with_binary_op("__ror__", other)
+
+    def __xor__(self, other: Any) -> Expression:
+        return self._with_binary_op("__xor__", other)
+
+    def __rxor__(self, other: Any) -> Expression:
+        return self._with_binary_op("__rxor__", other)
+
+    def __invert__(self) -> Expression:
+        return Expression(lambda df: ~self(df), f"(~{self._repr_str})")
+
     def __array_ufunc__(
         self, ufunc: Callable[..., Any], method: str, *inputs: Any, **kwargs: Any
     ) -> Expression:
diff --git a/pandas/core/ops/array_ops.py b/pandas/core/ops/array_ops.py
index 7b21772b443f6..d179a280d3e55 100644
--- a/pandas/core/ops/array_ops.py
+++ b/pandas/core/ops/array_ops.py
@@ -113,6 +113,9 @@ def fill_binop(left, right, fill_value):
 
 
 def comp_method_OBJECT_ARRAY(op, x, y):
+    from pandas._libs import missing as libmissing
+    from pandas.core.arrays import BooleanArray
+
     if isinstance(y, list):
         # e.g. test_tuple_categories
         y = construct_1d_object_array_from_listlike(y)
@@ -129,7 +132,31 @@ def comp_method_OBJECT_ARRAY(op, x, y):
         result = libops.vec_compare(x.ravel(), y.ravel(), op)
     else:
         result = libops.scalar_compare(x.ravel(), y, op)
-    return result.reshape(x.shape)
+    result = result.reshape(x.shape)
+
+    # GH#63328: Check if there are pd.NA values in the input and return
+    # BooleanArray to properly propagate NA in comparisons
+    x_has_na = any(val is libmissing.NA for val in x.ravel())
+    y_has_na = (
+        is_scalar(y) and y is libmissing.NA
+    ) or (
+        isinstance(y, np.ndarray)
+        and any(val is libmissing.NA for val in y.ravel())
+    )
+
+    if x_has_na or y_has_na:
+        # Create a mask for NA values
+        mask = np.array([val is libmissing.NA for val in x.ravel()], dtype=bool)
+        if isinstance(y, np.ndarray):
+            mask = mask | np.array(
+                [val is libmissing.NA for val in y.ravel()], dtype=bool
+            )
+        elif y is libmissing.NA:
+            mask = np.ones(x.shape, dtype=bool)
+        mask = mask.reshape(x.shape)
+        return BooleanArray(result, mask, copy=False)
+
+    return result
 
 
 def _masked_arith_op(x: np.ndarray, y, op) -> np.ndarray: