On the derivation of confidence intervals for hypothesis testing I understand that a 95% confidence interval should contain the value of the population parameter under the null hypothesis 95% of the time. That is, for $H_0: \theta = \theta_0$:
$$
\begin{align*}
&\Pr(|t| \leq t_{crit}) = 0.95 \text{ under } H_0.\\
\Rightarrow &\Pr\left(\frac{\hat{\theta} - \theta_0}{se(\hat{\theta})} \leq t_{crit}\right) = 0.95 \text{ under } H_0.\\
\Rightarrow &\Pr\left(\hat{\theta} - t_{crit} se(\hat{\theta}) \leq \theta_0 \leq \hat{\theta} + t_{crit} se(\hat{\theta})\right) = 0.95 \text{ under } H_0.
\end{align*}
$$
That is an interval centered around $\theta_0$. Another way to write it down is:
$$
CI = \left[\hat{\theta} - t_{crit} se(\hat{\theta}), \hat{\theta} + t_{crit} se(\hat{\theta})\right]
$$
We can then reject $H_0$ (at a 5% significance level) if $\theta_0$ is not an element of $CI$. My questions is how can that be the case if by construction, we are centering $CI$ precisely around $\theta_0$.

Demetri Pananos' post already does the job, IMO.
Nevertheless, it's tempting to leave behind an added clarification as to what CIs do, as is done in Casella, Berger.
While working with interval estimators, for the matter random intervals $[L(\mathbf X),U (\mathbf X) ], ~L(\mathbf X)\leq U (\mathbf X), ~\mathbf X\in\mathcal X, $ it must be borne in mind that probability statements pertinent to the interval i.e. $\mathbb P_\theta(\theta\in [L(\mathbf X),U (\mathbf X) ]) $ "refer to $\mathbf X,$ not $\theta. $" More elaborately, here the interval is the random quantity and not the parameter: the interval has the probability to cover the parameter, and it is definitely not the case that the probability is about whether the parameter is inside the interval.
For illustration, let's talk about the very familiar CI of $\mu$ of $\mathcal N(\mu, \sigma^2),~\sigma $ being known. The $(1-\alpha) $ CI for $\mu$ is
$$\left[\bar X -\frac{z_{\alpha/2}\sigma}{\sqrt n}, ~\bar X+ \frac{z_{\alpha/2}\sigma}{\sqrt n}  \right]\tag 1.$$
There is not $100(1-\alpha) \%$ chance that $\mu$ is in the realised value of the interval $\left[\bar x -\frac{z_{\alpha/2}\sigma}{\sqrt n}, ~\bar x+ \frac{z_{\alpha/2}\sigma}{\sqrt n}  \right]: \mu$ is fixed; $\mu$ would either lie in the realised interval with probability $1$ or $0$ otherwise. The correct interpretation, which actually provides a better insight, is that, if repeatedly sampled from $\mathcal N(\mu, \sigma^2) $ (not necessarily samples of same size), then $100(1-\alpha) \%$ of the sampled points of the random interval $(1) $ cover the true $\mu.$

Reference:
Statistical Inference, George Casella, Roger L. Berger,
Duxbury, 2002.
