Using multiple imputation to estimate counterfactual means and a treatment effect is a valid way to do so and relies on the same assumptions that other causal inference methods that rely on covariate adjustment do. This strategy and the more common but similar strategy of g-computation are described in Westreich et al. (2015). It works as follows:
Create two new variables in your dataset, $Y^0$ and $Y^1$. For all units with $T=0$ (where $T$ is the treatment), set $Y^0=Y$ (where $Y$ is the observed outcome), and leave the rest as missing. For all units with $T=1$, set $Y^1=Y$, and leave the rest as missing. Remove $Y$ from your dataset. You can now perform multiple imputation as usual (e.g., using MICE), and the causal effect in each dataset is estimated as the mean of the difference between the two variables (i.e., $Y^1$ and $Y^0$). You can fit a regression of Y1 - Y0 ~ 1 in each imputed dataset, and the coefficient on the intercept will be the treatment effect estimate. You can then use Rubin's rules to pool the estimated standard errors and effect estimates across datasets to arrive at a single estimate and its standard error. It may also be a good idea to bootstrap instead of relying on the pooling formula for the standard error. Bartlett and Hughes (2020) describe bootstrapping approaches appropriate for use with multiple imputation.
This procedure is not a common way to estimate causal effects. There are two more common methods: regression imputation and matching imputation.
Regression imputation, also known as g-computation and the g-formula, is described in Snowden et al. (2011) and involves fitting a model for the relationship between the covariates, the treatment, and the outcome. From this model, predicted values of the outcome are generated for all units setting $T$ to 1 and then setting $T$ to 0, which yields estimated potential outcomes under each treatment for each unit. Again, you can impute the treatment effect estimate as the mean of the difference between the two potential outcomes across units. When using a parametric model for the outcome (e.g., linear or logistic regression), there are analytic approximations to the standard error of the treatment effect; this is implemented in Stata using teffects ra. Otherwise (i.e., when using a machine learning model for the outcome), valid inference is only possible when using a doubly-robust estimator (e.g., AIPW or TMLE), which additionally involves a model for the treatment given the covariates, and bootstrapping; see Daniel (2014) for an introduction.
Matching imputation is described in Abadie and Imbens (2006) and involves finding a control unit that is close to each treated unit and using the control unit's outcome as the potential outcome under control for that treated unit. You can do the same by finding a treated unit that is close to each control and using the treated unit's outcome as the potential outcome under treatment for that control unit. There are many ways to define "closeness" that have different properties and many ways to customize the matching, such as by using multiple matches and taking the average or restricting the allowable distance between each unit and its match(es). Matching imputation doesn't require a model for the outcome, which can be appealing, though its performance can be improved when a model for the outcome is incorporated as described in Abadie and Imbens (2011). This method is implemented in teffects nnmatch in Stata. Matching can also be used for causal effect estimation without performing imputation; this philosophy of matching as "nonparametric preprocessing" is described by Ho et al. (2007) and uses matching to mimic the conditions of a randomized trial by equating the distributions of the covariates in the treated and control units.
There are other methods of estimating causal effects using covariate adjustment, but those are the ones based on imputation. You might also look into propensity score weighting, another popular approach that does not involve imputation (but can be combined with g-computation).
