Why is the median of an even number of samples the arithmetic mean? When calculating the median for an even number of samples, we sort the samples and take the arithmetic mean of the two in the middle.
This seems quite arbitrary - it assumes the best way to interpolate between the two samples is a linear interpolation.
Why not use the geometric or some other mean? Wikipedia lists quite a lot of options:

I guess one argument is that a choice needs to be made and any choice would add some bias. Is there a more fundamental reason?
 A: When this happens, I see two options.

*

*Report the median as an interval, which gives us infinitely-many medians instead of just one. (This is not even a confidence interval or credible interval for a point estimate.)


*Figure out a way to pick one value.
Picking the arithmetic mean of the middle two points is just one way of doing the estimation. Advantages include the ease with which is is computed and explained. Disadvantages could include bias when there is a skewed distribution.
There are always many ways of estimating a quantity. In fact, the quantile function in R software has at least nine options for computing quantiles like the median. If you have reason to believe that an arithmetic mean of the middle two values has inferior properties to some other form of estimation, you are free to argue why your alternative is better. Plenty of people have seen other estimators that are not optimal for them and have proposed alternatives (e.g., James-Stein estimator).
EDIT
The following simulation shows the claim by John Madden in the comments that all values in the interval minimize the absolute deviation function whose minimizer is one way to define median.
set.seed(2023)
y <- c(1, 2, 3, 7, 8, 9)
candidate_medians <- seq(2, 8, 0.01)
median_loss <- rep(NA, length(candidate_medians))
for (i in 1:length(candidate_medians)){
  
  median_loss[i] <- mean(abs(y - candidate_medians[i]))
  
}
plot(candidate_medians, median_loss)


All values $M$ between the middle data values, $3$ and $7$ (so $M\in[3,7]$), give equal and minimizing values of $L(y, M) = \frac{1}{n}\sum_{i=1}^n\big\vert y_i-M\big\vert$.
A: One of the first and best uses of the median is to identify the location (central value) of a dataset or distribution in a robust way.  From this perspective it doesn't matter how you choose among alternative values when there is more than one median: they are all equally good.
However, the very meaning of "central value" is questionable when a distribution is not symmetric.  (See Why is median age a better statistic than mean age? for some discussion.) But it does have a meaning when the distribution is symmetric near its middle.  (The sense of "distribution" includes any dataset, which we may consider to be equivalent to its empirical distribution.)
This means there is some number $\epsilon \gt 0$ for which the distribution is symmetric, or at least very nearly so, for all values within a distance $\epsilon$ of its median.  Since "distance" if measured on the scale of the data would be data-dependent, let's measure it in terms of probability.  That is, let us say a distribution $F$ (the cumulative function) is "symmetric near a median $\tilde \mu$" when there exists $\epsilon \gt 0$ for which $|F(\tilde\mu+\delta)-1/2|\le \epsilon$ implies $F(\tilde\mu+\delta) + F(\tilde\mu-\delta)=1.$  See https://stats.stackexchange.com/a/29010/919 for a more precise and general approach using a similar formula.  The idea is that the shape of the middle part of the distribution comprising a total probability (or proportion) of $2\epsilon$ is symmetric, but using values $\epsilon \lt 1/2$ permits asymmetry to appear in the tails of the distribution.
For discrete distributions (such as empirical distributions) perfect symmetry will rarely be the case due to the jumps in $F$ of size $1/n$ for datasets of size $n.$  We should therefore be willing to relax this requirement to "approximately symmetric in its center" when we have chosen a value $\lambda \ge 0$ and a median $\tilde\mu$ for which
$$|F(\tilde\mu+\delta)-1/2|\le \epsilon \text{ implies }|F(\tilde\mu+\delta) + F(\tilde\mu-\delta)-1| \lt \frac{\lambda}{n}$$
$\lambda$ measures the degree of asymmetry and a value of $\lambda = 1$ (or even a little larger) ought to be acceptable.
One basic technique of data analysis is re-expression (also known as transformation).  Depending on the application and meaning of the data, we seek a "natural" or "simple" function $f,$ such as a power (or Box-Cox transformation), which makes the re-expressed data approximately symmetric in its center.  If this can be achieved, it permits simplified descriptions and attendant insights, because we can summarize the distribution by giving a single central value that accounts for much of what the data are doing.  After that we can focus on describing how its tails might depart from the central part of the distribution: how long or heavy they appear and their degree of asymmetry.
Once we have found such a transformation (it needn't be unique) and applied it to re-express the data, the arithmetic mean of the extreme possible values of the median is a meaningful and appropriate value for the center of a distribution, because it is invariant under reversal of the data -- which, due to the approximate central symmetry, does not appreciably alter the shape of the central part of the distribution.  Thus, any proposed median exceeding the arithmetic mean is just as valid as a proposed median falling equally far below the arithmetic mean.  The middle is the intuitively obvious choice.
John Tukey exploits a generalization of this idea to develop robust ways to identify and measure the asymmetry of any positive distribution and select a Box-Cox parameter that will make it approximately symmetric in its center.  Writing $F_{-}$ for the distribution of the reversed data (that is, $F_{-}(x)$ is the proportion of the data equal to or greater than $x$), you study how the mid-quantiles $\mu_q(F)=(F^{-1}(q) + F_{-}^{-1}(q))/2$ vary with $q,$ for a suitably chosen series of $q$ that probe the tails: usually $q = 1/2, 1/4, 1/8, 1/16, \cdots.$  These mid-quantiles will not vary for as long as the central $1-2q$ portion of the distribution is symmetric.  The arithmetic mean of the two middle values in an even-size dataset is precisely $\mu_{1/2}(F).$  For the details, see Tukey's book EDA (Addison-Wesley 1977).  My post at https://stats.stackexchange.com/a/582120/919 provides a worked example.
