We can construct the regression:
$$\begin{bmatrix} y_1 \\ y_2 \end{bmatrix} =
\begin{bmatrix} x_1 & 0 \\ x_2 & x_2 \end{bmatrix}
\begin{bmatrix} \beta_1 \\ \Delta \end{bmatrix} +
\begin{bmatrix} e_1 \\ e_2 \end{bmatrix}$$
When estimating the parameters, we should keep in mind that $\sigma^2_1$ need not equal $\sigma^2_2$, which introduces a simple, well-structured, heteroskedasticity into the estimation, which can be addressed via weighted least squares.
The null hypothesis is $H_0: \Delta = 0$, and the alternative is, evidently, $H_A: \Delta \neq 0$. In the case of Gaussian errors, the obvious test is an $F$-test, which would be exact in finite samples if it were not for the likely mild) disruption caused by having to estimate two variance terms instead of one.
To see how much of a disruption the heteroskedasticity causes to the distribution of the $F$-statistic, we construct an example with $n_1 = n_2 = 100$ and four regressors in both $x_1$ and $x_2$. We set $\sigma^2_2 = 4\sigma^2_1$, and estimate the regression using iteratively reweighted least squares. We then calculate the p-value of the $F$-test, which, under $H_0$, should be distributed according to a Uniform distribution. Repeating the entire process 10,000 times allows us to test the distribution of the 10,000 p-values against the Uniform distribution, in this case using a Kolmogorov-Smirnov test:
p_values <- rep(0,10000)
Zeroes <- matrix(0, nrow=100, ncol=4)
for (j in 1:length(p_values)) {
X1 <- matrix(rnorm(400), nrow=100)
X2 <- matrix(rnorm(400), nrow=100)
y1 <- rowSums(X1) + rnorm(100)
y2 <- rowSums(X2) + 2*rnorm(100)
y <- c(y1, y2)
X <- rbind(cbind(X1, Zeroes),
cbind(X2, X2))
# Iteratively reweighted least squares
repeat {
ehat <- residuals(lm(y~X))
var_ratio <- var(ehat[101:200]) / var(ehat[1:100])
y[101:200] <- y[101:200] / sqrt(var_ratio)
X[101:200,] <- X[101:200,] / sqrt(var_ratio)
if (abs(var_ratio-1) < 0.00001) break
}
m1 <- lm(y~X)
m2 <- lm(y~X[1:200,1:4])
p_values[j] <-anova(m1, m2)$`Pr(>F)`[2]
}
ks.test(p_values, "punif")
One-sample Kolmogorov-Smirnov test
data: p_values
D = 0.008737, p-value = 0.43
alternative hypothesis: two-sided
which indicates, at least in this case, that the $F$-statistic's distribution is pretty close to the theoretical distribution, even for a not-terribly-large sample size.
In the case of non-Gaussian errors, it seems not unlikely that the $F$-test will break down, due to its known sensitivity to exactly this situation. In that case, an alternative would be to construct a permutation or bootstrap test (https://en.wikipedia.org/wiki/Resampling_(statistics)) based on the $F$-statistic, but I should point out that if you were to do so, there would likely be no particular advantage to sticking with the $F$-statistic as the test statistic of choice. An asymptotically equivalent permutation test can be applied if $n_1$ and $n_2$ together are too large for an exact permutation test; the asymptotics are such that the proper significance level is obtained (asymptotically) as long as exchangeability assumptions are met (and of course the null hypothesis is true.)
You can also rely on the asymptotic distribution of the $F$-statistic under non-Gaussianity, as outlined here: Does the F-test for multivariable regression work with non-normal residuals but large sample size?. However, IIRC the power of the $F$-test can be severely affected even if its significance level is approximately correct, hence my recommendation for the use of a permutation-based test.