The problem can be solved if the distribution $f(\cdot)$ is in a Levy stable distribution family. In your concrete example, since the $d\dim$ normal distribution $N(\mu,I_d)$ is "additive", the exact distribution of the $k$-sum $X_{i_1}+\cdots+X_{i_k}$ is $N(k\mu,kI_d)$. The choice of these $k$ random variables $X_{i_1},\cdots,X_{i_k}$ does not affect the mean nor covariance since the $X_i$'s are i.i.d. So it suffices to derive the distribution of $\|Y\|^2$ for a $d$-dimensional $Y\sim N(k\mu,kI_d)$, and there are $\left(\begin{array}{c}
n\\k\end{array}\right)$ possible combinations of one realization value of $\|Y\|^2$ if $Y=X_{i_1}+\cdots+X_{i_k}$ is the sum of an i.i.d. sample of size $k$ therefore the probability $Pr(max\|X_{i_1}+\cdots+X_{i_k}\|^2=y)$ will be $\left(\begin{array}{c}
n\\k\end{array}\right)\cdot f_{\|Y\|^2}(y)\cdot F_{\|Y'\|^2}(y)$. The maximum is taken over all possible $1 \leq i_{1} < i_{2} < \ldots < i_{k} \leq n$, same below. The $Y'=\sum_{j\notin\{i_1,\cdots,i_k\} }X_j$ follows a similar derivation but $F$ is the c.d.f.
To calculate the latter $f_{\|Y\|^2}(y)$, we need to notice that $\frac{Y}{\sqrt{k}}\sim N_d(\sqrt{k}\mu,I_d)$, we can compute its $L^2$ norm $\frac{1}{k}\sum_iY^2_i=:\|\frac{Y}{\sqrt{k}}\|^2\sim\chi^2(d)$. Thus the p.d.f. $f_{\|Y\|^2}(y)=f_{\chi^2(d)}(ky)$
To draw a conclusion, $Pr(max\|X_{i_1}+\cdots+X_{i_k}\|^2=y)=\left(\begin{array}{c}
n\\k\end{array}\right)\cdot f_{\chi^2(d)}(ky)\cdot F_{\chi^2(d)}((n-k)y)$ in case of $d\dim$ multivariate normal as you indicated in OP. For other families of $f(\cdot )$ closed under addition, we can proceed the same way, and the resulting density may not be of closed form if we do not know one under square transformation (in contrast that we know $\chi^2(d)$ in above.).
