After revisiting my question, I think I have managed to find a proof that there are no other examples of lattices with the requested property. I'm posting it as a self-answer in case someone is interested. Here is the outline:

 * The case of even dimension $\mathrm{dim}\: \Lambda \ge 4$ was already described$^\dagger$ in the "Attempt" section of my question; the most important part (bounding the generalized Bernoulli numbers) is done in Section 3 of [this 2011 paper](https://arxiv.org/abs/1110.6430) by M. Johnson which I cited in the question.

 * We can also adapt the strategy in that paper to the odd-dimensional case. As in the even case, we start by observing that the theta series must be a certain Eisenstein series of half-integer weight $k$ and prescribed level $4N$ and character $\chi$, normalized such that the zeroth coefficient is $1$. One needs an expression for the Fourier coefficients of these series, which I found for $\mathrm{dim}\: \Lambda = 2k \ge 5$ in page 17 of [this 2015 thesis](https://escholarship.org/content/qt98k5g03d/qt98k5g03d.pdf) by M. Owen. We are only interested in the coefficient of $q^1$, which simplifies to

  $$c_1 = \left(-2\pi i/N\right)^k \:\Gamma(k)^{-1} A(1) X(1),$$

  where

  $$A(1) = \sum_{r=1}^{4N} \epsilon_r^{2k} \left(\frac{4N}{r}\right) \chi(r) e^{2\pi ir/N}, \quad X(1)= \prod_{l\nmid 4N}(1+\chi(l) \epsilon_l^{2k-1} l^{1/2-k})$$

  (as in Owen's paper, here $l$ is prime, $\left(\frac{4N}{r}\right)$ is a Kronecker symbol, and $\epsilon_n$ denotes the principal branch of $\sqrt{\left(\frac{-1}{n}\right)}$. Note also that I use $4N$ where the author uses $N$, since for half-integer weight forms the level is always a multiple of 4). This $c_1$ is an analogue of the generalized Bernoulli numbers for half-integral weight.

  The $A(1)$ factor is a sum of roots of unity whose modulus is obviously bounded by $4N$, and by the same argument as in page 9 of Johnson's paper, $|X(1)|$ is bounded by $\prod_{l}(1+l^{1/2-k})=\frac{\zeta(k-1/2)}{\zeta(2k-1)}$. Putting everything together, we see that the bound on $|c_1|$ will be lower than $2$ for high enough $N$ and $k$ (concretely we have either $N=1$ and $k \le 9/2$ or $1 < N \le 3$ and $k=5/2$), so recalling that $c_1$ must be a nonzero even integer, we can restrict our analysis to these cases only. All of them are then ruled out by manual computation of the coefficients.

 * The case $\mathrm{dim}\: \Lambda = 2$ is easily dealt with by using the relationship between 2D integer lattices and rings of integers of imaginary quadratic number fields, and applying the class number formula.

 * This only leaves the case $\mathrm{dim}\: \Lambda = 3$. It is known that there are no isospectral lattices in dimension lower than $4$ (see e.g. the end of Section 2 [here](https://www.maths.ed.ac.uk/~v1ranick/papers/conwaysens.pdf)), so in this case it suffices to appeal to the Siegel-Weil formula and check the single-class examples in dimension $3$ with the LMFDB online tool, which I already did (this could also be an alternative way to treat the case $\mathrm{dim}\: \Lambda = 2$). This completes the proof.

----------

**Remark:** I admit that the proof above is somewhat "ugly", since it does not give any insight into the form of the classification. A more satisfying proof would explain why all of these lattices have the multiplicative structure of an order in a division algebra (and thus occur in dimensions 1, 2, 4, 8) except for the sporadic case $(E_8 \times E_8, D_{16}^+)$.

There is an observation to be made regarding that last case. For any lattice $\Lambda$ it is possible to define [Siegel theta series](https://en.wikipedia.org/wiki/Siegel_theta_series) of genus $g$ generalizing the usual theta series, which essentially count the number of $g$-dimensional sublattices of $\Lambda$. Since there is a bijective correspondence between lattices and their theta series in genus $g\ge \mathrm{dim} \: \Lambda$, the issue cannot arise that two lattices have the same theta series; in fact, the sporadic example of $E_8 \times E_8 \leftrightarrow D_{16}^+$ disappears already in genus 4 because of the [Schottky form](https://en.wikipedia.org/wiki/Schottky_form), and the corresponding Siegel-Weil formula for $g\ge 4$ takes the form $\mathrm{E}_8^{(g)} = \frac{405}{691} \Theta_{E_8 \times E_8}^{(g)} + \frac{286}{691} \Theta_{D_{16}^+}^{(g)}$ with $\Theta_{E_8 \times E_8}^{(g)} \neq \Theta_{D_{16}^+}^{(g)}$. In contrast, the octonionic "$\mathrm{E} = \Theta$" identity between the Siegel theta series of the $E_8$ lattice and the Siegel Eisenstein series of weight 4 does hold for all genera, and I believe the same happens for the other examples. This could perhaps be a starting point to explain why the sporadic case doesn't have a multiplicative structure.

----------

$\dagger$- The extra possible example of dimension $4$ and level $42$ that I mentioned does not correspond to any lattice; perhaps the easiest way to see this is to inspect the coefficients of the candidate theta series, which starts as

$$1+2q+2q^2+2q^3+2q^4+12q^5+2q^6+2q^7+2q^8+2q^9+12q^{10}+\ldots$$

(the general coefficient is twice the sum of the divisors of $n$ that don't divide $42$).

Since the coefficient of $q^n$ counts the number of vectors of norm $2n$, such a lattice would have a vector $\mathbf{a}$ of squared norm $2$ and another, linearly independent vector $\mathbf{b}$ of squared norm $10$. We can also see that all vectors of squared norm lower than $10$, as well as those with squared norms strictly between $10$ and $20$, must be integer multiples of $\mathbf{a}$, since there are only two vectors of each norm in these intervals.

The squared norm of the sum $\mathbf{a}+\mathbf{b}$ is a positive even integer bounded by $(\sqrt{2}+\sqrt{10})^2 \approx 20.944$; by the above observations, it can only be $20$. But then by the parallelogram law, $|\mathbf{a}-\mathbf{b}|^2 = 2|\mathbf{a}|^2+2|\mathbf{b}|^2-|\mathbf{a}+\mathbf{b}|^2 = 4+20-20=4$, which is a contradiction since $\mathbf{a}-\mathbf{b} \not\propto \mathbf{a}$.