how to perform double ML with binary data (either in the treatment or in the outcome)? I have grown interested in double Machine Learning (ML) for causal inference because it answers an intuitive question: if the relationship between a variable $X$ (the treatment) and a variable $Y$ (the outcome) is confounded by a third variable $Z$ (the confounder), can't we extract some unconfounded signal to get an unbiased estimate of $X$ on $Y$?

That's what double ML does! By extracting the residuals $\epsilon_{X}$ and $\epsilon_{Y}$ (from the regressions of $X$ and $Y$ on $Z$), we "purify" the effect of $Z$ on $X$ and $Y$. We can then estimate the effect of $X$ on $Y$ from those residuals which are, by construction, unconfounded.
However, as I wanted to apply this at work, I rapidly stumbled upon another difficulty: how can we apply this approach when the treatment $X$ or the outcome $Y$ is a binary variable?
In this book, it is recommended to not use double ML in a binary outcome context (but I even fail to see how it can implemented) : "If the treatment or the outcome is binary, one might think it is better to replace the machine learning regression models for their classification versions. However, this does not work. The theory of orthogonalization only functions under regression models, similarly with what we’ve seen a long time ago when talking about Instrumental Variables. To be honest, it is not that the model will fail miserably if you replace regression by classification, but I would advise against it."
Is there any workaround solution in this case?
 A: There are no issues with using the double ML strategy provided that your classifier provides probabilistic outputs, i.e. it estimates the conditional response function. However in finite samples, it is possible the strategy produces results outside the natural bounds of the problem. For example, the average treatment effect with a binary outcome must fall within [-1,1].
This happens because double ML supplies a bias correction term to a plug-in estimate of your target estimand, and this correction term does not necessarily respect the constraints of the problem in a finite sample. However, it practically almost always does, so I would not be afraid to use this estimation methodology - plus, the theory still works out in this case.
An alternative approach you could try is TMLE. The two methodologies are quite similar in their flexible use of ML estimators for nuisance functions and the ability to use double ML'S idea of cross-fitting (called CV-TMLE in this literature).
However, TMLE has an advantage in this case in that it does infact respect problem constraints in finite samples. It achieves this by providing a correction in an earlier step. We still use a plug-in estimator for our target estimand, but now we provide the correction to the quantity being plugged-in. The relevant quantity being plugged-in in this case would be the log-odds of the outcome under treatment or control. Since we provide the correction on the log-odds scale, when we convert back to the probability scale when we plug-in, we respect the bounds of the problem.
