You don't say anything about the populations from which your data were randomly sampled. If sample sizes are at least moderately large and the data do not show
signs of severe skewness (e.g., far outliers on one one side of center), then
a two sample t confidence interval is probably the simplest and best choice.
If data are from continuous non-normal populations, then a two sample Wilcoxon (rank sum) procedure may provide a useful confidence interval (CI).
For example, suppose we have samples of sizes $n_1 = 30$ and $n_2 = 20$ from respective distributions $\textsf{Gamma}(\text{shape}=5, \text{rate}=.2)$ and
$\textsf{Gamma}(6, .4).$
set.seed(2019); n1=30; n2=20
x1 = rgamma(n1, 5, .2); x2 = rgamma(n2, 6, .4)
summary(x1); sd(x1)
Min. 1st Qu. Median Mean 3rd Qu. Max.
7.653 16.314 24.456 27.025 33.651 59.145
[1] 13.13563 # sample SD
summary(x2); sd(x2)
Min. 1st Qu. Median Mean 3rd Qu. Max.
5.847 9.563 14.717 14.795 19.022 26.130
[1] 6.076901
In in each sample, the sample mean and median are near each other.
Boxplots show no outliers. Also, notches in the boxes are nonparametric CIs; non-overlapping CIs suggest that the two population
medians are different.
boxplot(x1, x2, notch=T, col="skyblue2")

Especially knowing the populations, one might raise theoretical objections
to either the 95% Welch t CI $(6.67, 17.79)$ or the 95% Wilcoxon CI $(4.79, 17.11)$. But we know less just looking at descriptions of the sample, and there are no strong warnings in those descriptions against using either kind of CI. Results from relevant R procedures are shown below. Neither CI includes $0,$ so one can conclude that the difference
$\delta = \mu_1 - \mu_2$ is significantly different from $0.$
t.test(x1, x2)$conf.int
[1] 6.673883 17.786275
attr(,"conf.level")
[1] 0.95
wilcox.test(x1, x2, conf.int=T)$conf.int
[1] 4.785996 17.109673
attr(,"conf.level")
[1] 0.95
However, you specifically asked about making a bootstrap confidence
interval for the difference in population means. Your interest in making a bootstrap CI may be that you want to have the experience doing one, or that doing a bootstrap is a class assignment. There are very many
styles of bootstrap CIs; I will illustrate one of the simpler styles.
You seek a CI for $\delta = \mu_1 - \mu_2,$ which can be estimated
by $D = \bar X_1 - \bar X_2.$ If you knew the distribution of $D,$ then you could use it to find bounds $L$ and $U$ such that
$$0.95 = P(L \le D - \delta \le U) = P(-U+D \le \delta \le -L+D),$$
so that a 95% CI for $\delta$ would be of the form $(D-U, D-L).$
Recall that our observed value of $D$ is $D_{obs}= \bar X_1- \bar X_2 = 12.23,$
which we use in the bootstrap procedure.
d.obs = mean(x1) - mean(x2); d.obs
[1] 12.23008
Not knowing the distribution of $D$ we use bootstrapping to find
useful approximations $L^*$ of $L$ and $U^*$ of $U.$ In order to do this we take a large number $B = 2000$ of samples of size 30 with replacement
form x1
and then find their means $A_1^*.$ Similarly, we 're-sample' from x2
to find means $A_2^*.$ From these, we find a vector of 2000 quantities $D^* = A_1^* - A_2^*.$ We would like to have a 'bootstrap distribution' of $D - \delta.$ Because $\delta$ is unknown, we use
$D_{obs}$ as a proxy.
Then we find $L* \approx L$ and $U^* \approx U$
by taking quantiles .025 and .975 of $D^* - D_{obs}.$ From there, the 95% bootstrap Ci for $\delta$ is of the form
$(D_{obs}-U^*, D_{obs}-L^*),$ where $D_{obs}$ returns to its original role.
In the R code below, we use .b
instead of $*$ to indicate quantities arising from re-sampling. The 95% bootstrap CI for $\delta$ is
$(6.99, 17.29).$ This is called a nonparametric CI because the
bootstrap procedure did not involve the assumption that data are normal.
set.seed(1234); B = 2000
D.b = replicate(B, mean(sample(x1,n1,rep=T)) -
mean(sample(x2,n2,rep=T)) - d.obs)
d.obs - quantile(D.b, c(.975,.025))
97.5% 2.5%
6.991987 17.289737
Addendum (per revision of question and comment): I agree with @Dave, that the t.test
should be OK for such a large sample--even if markedly skewed.
I do not
have much experience bootstrapping huge samples,
but the bootstrap works well in the example I
show below, comparing large exponential samples with respective population means 100 and 102.
The 95% t CI is $(1.016. 2.343)$ and the 95%
nonparametric bootstrap CI is $(1.026, 2.344).$
set.seed(807)
x1 = rexp(10^5, 1/102); x2=rexp(10^6, 1/100)
t.test(x1, x2)$conf.int
[1] 1.016061 2.343464
attr(,"conf.level")
[1] 0.95
set.seed(1234); B = 2000
n1 = length(x1); n2 = length(x2)
d.obs = mean(x1) - mean(x2)
D.b = replicate(B, mean(sample(x1,n1,rep=T)) -
mean(sample(x2,n2,rep=T))-d.obs)
d.obs-quantile(D.b, c(.975,.025))
97.5% 2.5%
1.02571 2.34406