Simulating MADMAX in 3D: Requirements for Dielectric Axion Haloscopes

We present 3D calculations for dielectric haloscopes such as the currently envisioned MADMAX experiment. For ideal systems with perfectly flat, parallel and isotropic dielectric disks of finite diameter, we find that a geometrical form factor reduces the emitted power by up to $30\,\%$ compared to earlier 1D calculations. We derive the emitted beam shape, which is important for antenna design. We show that realistic dark matter axion velocities of $10^{-3} c$ and inhomogeneities of the external magnetic field at the scale of $10\,\%$ have negligible impact on the sensitivity of MADMAX. We investigate design requirements for which the emitted power changes by less than $20\,\%$ for a benchmark boost factor with a bandwidth of $50\,{\rm MHz}$ at $22\,{\rm GHz}$, corresponding to an axion mass of $90\,\mu{\rm eV}$. We find that the maximum allowed disk tilt is $100\,\mu{\rm m}$ divided by the disk diameter, the required disk planarity is $20\,\mu{\rm m}$ (min-to-max) or better, and the maximum allowed surface roughness is $100\,\mu{\rm m}$ (min-to-max). We show how using tiled dielectric disks glued together from multiple smaller patches can affect the beam shape and antenna coupling.

sioned MADMAX experiment. For ideal systems with perfectly flat, parallel and isotropic dielectric disks of finite diameter, we find that a geometrical form factor reduces the emitted power by up to 30 % compared to earlier 1D calculations. We derive the emitted beam shape, which is important for antenna design. We show that realistic dark matter axion velocities of 10 −3 c and inhomogeneities of the external magnetic field at the scale of 10 % have negligible impact on the sensitivity of MADMAX. We investigate design requirements for which the emitted power changes by less than 20 % for a benchmark boost factor with a bandwidth of 50 MHz at 22 GHz, corresponding to an axion mass of 90 µeV. We find that the maximum allowed disk tilt is 100 µm divided by the disk diameter, the required disk planarity is 20 µm (min-to-max) or better, and the maximum allowed surface roughness is 100 µm (min-to-max). We show how using tiled dielectric disks glued together from multiple smaller patches can affect the beam shape and antenna coupling.
In the presence of a strong external B-field, axions are converted into electromagnetic radiation at interfaces of media with different dielectric constants . The MADMAX experiment consists of a metallic mirror and many parallel dielectric disks in vacuum leading to electromagnetic radiation from each interface separating regions with different . Depending on the disk positions the radiation from different interfaces can interfere constructively and excite resonances between the dielectric disks, although with significantly lower quality factors as cavity experiments. The power boost factor β 2 describes the enhancement of the power emitted by the mirror together with the set of dielectric disks (booster) with respect to the radiation emitted by a perfect mirror of the same area and under the same B-field. Previous one dimensional (1D) calculations [22] showed that with 80 lanthanum aluminate (LaAlO 3 ) disks ( ≈ 24) a power boost factor of ≈ 5 × 10 4 can be achieved over a bandwidth of 50 MHz, leading to an emitted power of where A is the surface of the dielectric disks, B e the strength of the external magnetic field, ρ a the local cold dark matter density and |C aγ | a model-dependent coupling constant proportional to the axion-photon coupling g aγ as defined in [22], with typical values of |C aγ | ≈ 1.9 (KSVZ model [46,47]) or |C aγ | ≈ 0.7 (DFSZ model [48,49]). It is of central importance to understand the systematic uncertainties in the power boost factor β 2 . Previous work has relied on a 1D model for β 2 [22], while three dimensional (3D) effects have only been taken into account for smaller systems with up to one dielectric disk [50,51]. The work presented here extends these studies to systems with multiple dielectric disks as envisioned for MADMAX. We present simulations taking some of the most important realistic boundary conditions for an open booster (disks surrounded by free space) into account, i.e., first of all the fact that the disks are of finite size (ideal 3D booster), but also implications from a finite axion velocity, magnetic field inhomogeneities, mechanical tolerances, imprecise disk geometries, tilts and tiled disks (non-ideal booster). To this end we apply the finite element method by using the azimuthal symmetry of the booster (2D3D FEM), as well as the Recursive Fourier Propagation method, both introduced in [50]. In addition, we use the Mode Matching formalism briefly described in the next chapter. For a comparison showing their consistency see appendix A.1.
The paper in large parts is based on results from two PhD theses [52,53]. It is structured as follows: In section 2 we identify eigenmodes independently propagating inside the system, which form the basis for our description of the booster. Section 3 deals with finite-diameter but perfectly parallel and flat dielectric disks to which we refer as the ideal 3D case. Finally in section 4 we study non-ideal effects including effects from a finite axion velocity, B-field inhomogeneities, disk tilts and surface inaccuracies. We also discuss dielectric disks glued together from smaller uniform patches (tiled disks).

System Modes
At first order in the axion-photon coupling g aγ the axion-Maxwell equations can be written (using natural units with = c = 1 and the Lorentz-Heaviside convention α = e 2 /4π) as a wave equation for the electric field E using time-harmonic fields as [50] where ω = 2πν = m a is the angular frequency and is the permittivity. The permeability is assumed to be µ = 1. The axion field a on the right hand side acts as a source of electric fields, through its coupling constant g aγ and external magnetic field B e . For and B e constant over lengths much larger than the free photon wavelength λ = 2π/ω [50,54,55], a solution is given by the axion-induced field  Main panel: Spatial field distribution. Side and bottom panels: Fields on the x and y axis (at y = 0, x = 0, respectively) (blue solid lines) and fields corresponding to a Gaussian beam [58] with a beam waist radius of w 0 ≈ ø/3 (gray dashed lines), discussed later in section 3. Right: Same as main panel on the left, but for higher modes with = 0.
The axion-induced field E a has a discontinuity at a boundary between regions with different and hence does not solve eq. (2.1) anymore. The full solution is obtained by adding emitted electromagnetic radiation from the boundary compensating the discontinuity [22]. In all figures throughout this paper the electric fields are shown in units of E 0 ≡ max |E a | at a fixed instant of time if not stated otherwise.
For simplicity (and when not using FEM methods as e.g. in section 4.4), we will neglect free charges in the following by setting ∇ · E = 0. This sets the second term in eq. (2.1) to zero and the equation separates into three independent wave equations for each component of E, i.e., it is sufficient to consider each component as a scalar field (scalar diffraction theory). This approximation is valid for a dielectric haloscope with sufficiently homogeneous disks and has explicitly been confirmed for the ideal system discussed in the next section, as we show explicitly in appendix A.1. However, the calculations below can also be easily generalized by solving for the modes of the vectorized equation, see e.g. [56,57].
To begin with, consider a cylinder of dielectric material with radius R (diameter ø = 2R) surrounded by vacuum forming a dielectric waveguide. In the limit of large radius R λ and large dielectric constant 1, we obtain a model for one of the disks of the dielectric haloscope. In this limit the electric fields drop to zero at the outer boundary r = R of the disk. Explicitly, the solutions to the source-free scalar wave equations (i.e., eq. (2.1) with ∇·E = 0 and a = 0) are the eigenmodes which are illustrated in figure 1 and given by [56,57] 36 3 × 10 −5 . . . Table 1: Properties of the most important modes for a dielectric haloscope, as envisioned for the MADMAX prototype (disk diameter ø = 30 cm) and final scale experiment (disk diameter ø = 1 m). |η m | 2 denotes the coupling of the uniform axion field under a uniform magnetic field to the mode, k c,m its transverse momentum and δ d,m the diffraction loss parameter between the dielectric disks as defined in the text at 22 GHz.
with discrete radial mode indices m > 0 and azimuthal mode indices . J is the Bessel function of the first kind of order , and we take N m as a normalization factor such that |E m | 2 dA = 1. Here, k c is the transverse momentum, i.e., the momentum in the disk plane. These modes are orthogonal and complete in the sense that we can expand any field distribution inside of the disks into a set of these modes. Most importantly, they propagate independently along the z-direction within the disks as where e m are the coefficients for the mode expansion, k z,m is the propagation constant and k 0 = √ ω. In free space these eigenmodes of the dielectric disks in general do not propagate independently anymore, because they are no longer solutions of the scalar wave equation under the free space boundary conditions. Since they are orthogonal and complete, we still can expand fields at r < R into these modes, but during propagation they mix with each other, i.e., E(r, φ, z) = m, ,m , with the linear map P m m (z) between the modes. P can be calculated by using the scalar diffraction theory in free space discussed in [50]. P m m is the coefficient of the mode (m , ) when expanding the field obtained after propagating the mode (m, ) for a distance z in free space. One can generalize the 1D transfer matrix formalism for dielectric haloscopes in [22] by having left-moving and right-moving fields for each mode in each region and by directly including the mixing matrix P , see e.g. [52]. We refer to this kind of calculation as Mode Matching, because on an interface between two media with different dielectric constants the sum of modes describing the fields on one side needs to be matched with the respective sum of modes on the other side.
In order to see which modes are relevant for a dielectric haloscope, we have to consider their coupling to the axion-induced field E a . Table 1 summarizes the most important modes for the dielectric haloscope for a uniform external magnetic field and negligible axion velocity. The coefficients η m refer to the coupling of the mode to the axion-induced field E a , i.e., they are the coefficients of the modes when expanding E a on the disk surfaces into the modes. Explicitly, with a normalization factor N ( ) a such that |η m | 2 = 1. For our axion haloscope actually only the azimuthally symmetric ( = 0) lower modes with m = 1, 2, 3, 4 have a coupling stronger than 2% to the axion-induced electric field. All modes with = 0 do not couple due to symmetry, although imperfections may affect η, see section 4. When only considering these relevant modes with m = 1, 2, 3, 4; = 0 even for disks with diameter ø = 30 cm at 22 GHz the mixing between the modes, i.e., |P m m | for (m, ) = (m , ), is smaller than ≈ 8 × 10 −4 . So unless the system is tuned to be very resonant, the mixing can be neglected, i.e., P m m becomes diagonal and can be written as where δ d,m is a diffraction loss parameter and all modes propagate essentially independently. The parameter δ d,m can be suppressed by using disks with larger diameters, as expected.
If, on the other hand, the system is tuned to be very resonant for a specific mode, the difference in k z,m for the other modes will make them rapidly dephase, i.e., make all other modes irrelevant.

Ideal 3D Booster
We first consider an ideal but 3D booster with disks of finite extent, which are however still perfectly flat and parallel. We study two benchmark systems, tuned to an axion mass of m a ≈ 90 µeV (ν ≈ 22 GHz). The optimal boost factor bandwidth is given by a trade-off between disk readjustment time for tuning, and actual data taking time. The minimum bandwidth is further limited by losses. Here we consider a bandwidth of ≈ 50 MHz close to preliminary estimates of the optimal bandwidth maximizing scan speed for MADMAX [24]. We consider a booster with 20 lanthanum aluminate disks (assuming an isotropic dielectric constant of = 24) with a disk diameter of ø = 30 cm and thickness 1 mm as presently foreseen for the MADMAX prototype; in addition, we examine an 80 disk system with a disk diameter of ø = 1 m as envisioned in the final MADMAX setup [24,59]. All presented simulations assume free space surrounding for simplicity, see also [50,60]. This setup is expected to maximize diffraction losses. Detailed studies on the impact of using other different boundary conditions, e.g., conducting walls, will be discussed in future works. Figure 2 shows the power boost factor of such systems in terms of total emitted power (solid blue) and the power which can be coupled to an antenna receiving Gaussian beams as defined in [58] with beam waist radius w 0 ≈ ø/3 (dashed blue) compared to the 1D result (dashed gray). The double-peak or four-peak substructure, respectively, corresponds to different contributing resonances, for more details see [22]. Results from different numerical methods, i.e., 2D3D FEM, Recursive Fourier Propagation and Mode Matching are consistent up to percent level, which is negligible for the experiment's sensitivity to axion CDM and the axion-photon coupling |C aγ |. This confirms the validity of the scalar diffraction theory for Turning our attention to the results themselves, we first notice that the boost factor curve is shifted to higher frequencies compared to the 1D calculation. This is easily understood considering the phase evolution of the different modes along the booster. Due to the transverse momentum k c of the modes the phase changes slower along the z-direction compared to the 1D case according to eq. (2.4). Therefore, in order to have the same resonant behavior as in 1D one needs to "speed up" the phase evolution by going to slightly higher frequencies. For the lower modes (m, ) with small transverse momenta k c k z the frequency shift compared to the 1D calculation is where j ,m is the m-th zero of J , which roughly scales linearly with m. Since higher modes have higher transverse momenta, cf. table 1, the shift is more pronounced for higher modes.
As each mode propagates essentially independently through the system, no matter how the disk spacings in the system are tuned, for a fixed disk diameter the different modes always appear at the same frequency shifts relative to each other. The bandwidth above which higher modes start to become relevant is therefore ∆ν β ≈ ∆ν 2 0 − ∆ν 1 0 , which gives 55 MHz for the prototype booster and 5 MHz for the full-scale booster, consistent with figure 2.
Now considering the power emitted by the system, we see that in 3D the boost factor is reduced compared to the 1D calculations. Since all modes are orthogonal, the total power emitted is simply the sum of the power carried by each mode, as indicated by the stacked hatched regions in figure 2. In the benchmark case for the MADMAX prototype (left) we see that the second mode is already shifted by almost the full bandwidth of the boost factor itself and we essentially only get the power contributed by the first mode within a 50 MHz bandwidth. Since this mode couples to 69 % (independent of disk diameter) to the axion field, the boost factor is reduced by up to 30 % compared to the 1D case. This effect should be seen as a reduced coupling efficiency (form factor) of the system to the axion field and not as (diffraction) loss. Indeed, the diffraction loss of the first mode arising from the finite disk size in this case is smaller than δ d ≈ 10 −5 at this frequency (see table 1) which is negligible. This may not hold anymore when we consider the geometrical inaccuracies in section 4.
Lastly, we have to consider how to couple the power leaving the booster with an antenna into a receiver. The fundamental mode has a frequency-independent 97% matching ratio with a Gaussian beam [58] with a beam waist radius of w 0 ≈ ø/3, see also figure 1 (left). We therefore consider the coupling efficiencies to Gaussian beam antennas with this beam shape in this paper. Hence, in case only the fundamental mode contributes, we can achieve very good coupling efficiencies. In case the total power is also carried by higher modes, like in the 80 disk calculation in figure 2 (right), one can only receive significant power provided by the fundamental mode with the Gaussian antenna. This contribution is still 70 % of the total power due to the coupling of the axion field to the fundamental mode. However, small couplings of the higher modes to the Gaussian may interfere destructively when coupled to the antenna, further decreasing the received power. In principle it is possible to design an antenna which is matched to a more optimal combination of modes, as long as their relative phase stays roughly constant over the boost factor bandwidth -or in other words the total beam shape does not change drastically with frequency. For the initial stage of dielectric haloscopes this may already be a too elaborate approach. In summary, as long as the boost factor bandwidth is smaller than the difference between the frequency shifts of the first two modes, the optimal antenna is one that couples only to the fundamental mode. In particular, for the MADMAX prototype, designed for the frequency range from 18 to 24 GHz, an antenna system which couples to a Gaussian beam with beam waist radius of approximately 10 cm is close to optimal.
Not considered here, but crucial for a final experimental realization, might be possible reflections on the antenna, especially those of the higher modes, which after the reflection may couple and interfere destructively with the fundamental mode. For MADMAX such reflection effects have been already experimentally studied in [61]. There it was demonstrated on a 5-disk setup that adverse effects due to reflections may be significantly reduced by absorbing unwanted radiation in the vicinity of the antenna and calibrating out residual reflections using a dedicated model.

Non-Ideal Effects
A realistic system will always have inaccuracies, contrary to what was assumed in the previous section. Therefore, in the following we study the influence of axion velocity effects and inhomogeneities of the external magnetic field (causing changes to the axion-induced field E a ), as well as geometrical imperfections (tilts, planarity, surface roughness and tiling of the disks).

Axion Velocity
With non-zero axion velocity v a the axion field a and therefore also the axion-induced electric field E a acquire a spatial phase factor exp(−im a v a x) over the setup. A velocity along the booster axis causes phase differences between the disks and has been studied already in [62,63]. A transverse velocity v a, tilts the otherwise perpendicular angle of emission from the individual disks [64][65][66]. We can study the effect of this tilting by decomposing the axion-induced field into the above modes and observing how the coupling efficiencies η m change with transverse velocity. For the fundamental mode one finds analytically (4.1) This holds for small axion velocities, i.e., m a v a, R 1, which is applicable for the MADMAX boosters below about m a < 500 µeV. An exact result can be found in appendix B. Below m a = 100 µeV the effect on the full-scale MADMAX sensitivity is negligible. Although the effects may become more relevant for higher masses of m a = 100 − 400 µeV, still 90 % of power is left in the fundamental mode. In particular, one would have to take the average of the signal power over the CDM velocity distribution rather than just considering one velocity. Typical data-acquisition times for MADMAX before tuning to the next frequency band are expected to be at the order of a few days [24]. Since the earth rotates within the CDM 'wind', some of the velocity effects will average out and make the above reduction even milder.
The effect on the boost factor is explicitly demonstrated on the benchmark boost factor for the full-scale MADMAX booster in figure 3 (left) for an axion velocity exaggerated by up to one order of magnitude compared to realistic CDM velocities of v a ≈ 10 −3 c. We see that the received power (dashed lines) is degraded while the total power emitted by the haloscope (solid lines, almost on top of each other) remains almost unchanged, but is in the modes that do not couple to the antenna. No curve is shown for v a = 10 −3 c, since already v a = 2 × 10 −3 c does not significantly change the boost factor compared to the zero velocity case. For even higher axion velocities (not shown) nearly all power is contained in higher modes. Since they have higher k c the total power boost factor shifts to higher frequencies. Higher modes are more prone to diffraction losses and the inaccuracies of the setup described in the following sections. Therefore, also the total power emitted is reduced for higher velocities. For realistic CDM velocities of v a ≈ 10 −3 c, however, our benchmark boost factors are not changed significantly.
A finite axion velocity slightly tilts the emissions from individual disks. Therefore, the center of the beam emitted from the booster shifts away from the center of the disk as shown in figure 3 (right). This effect could in principle be used to build a velocity-sensitive haloscope after the discovery of the axion and to investigate and measure properties of the local dark matter halo, see for example [63,64].

Magnetic Field Inhomogeneity
Analogous to the velocity effects above, a transverse inhomogeneity of the magnetic field implies a corresponding inhomogeneity in the axion-induced field E a . Therefore, it changes the amount of power coupled into the different modes. For example a magnetic field proportional to the beam shape of the fundamental mode would cause the coupling efficiency of the first mode to be |η 10 | 2 = 100 %. Realization of such a magnet is, however, technically challenging and may increase magnet cost significantly. Here we consider a magnetic field amplitude with azimuthal and radial inhomogeneity, which is motivated by the symmetry of typically considered dipole magnets. We consider such a magnetic field parametrized by where B 0 is the magnetic field amplitude, h is the maximum relative scale of the inhomogeneity on the disk, R = ø/2 and k is a non-zero positive integer. For small h 1 one can show that the relative change in the coupling coefficients of the = 0 modes happens only at second order in h and is given by i.e., radial symmetric transverse inhomogeneities at the 10 % level leave the mode coupling coefficients unchanged well below the percent level. Therefore, such inhomogeneities have insignificant impact on sensitivity. This result has been confirmed with explicit numerical calculations using Recursive Fourier Propagation.

Geometrical Inaccuracies of the Dielectric Disks
Next we consider geometrical inaccuracies such as disk tilts, disk planarity and surface roughness. These mainly affect propagation of electromagnetic waves within the booster. If the distance between two interfaces varies as ∆z(r, φ) in the plane parallel to the disk surfaces (transverse thickness variation), the corresponding phase changes during propagation give rise to additional mode mixing as where we have left out the propagation and corresponding diffraction by a distance z in this formula for clarity (it is included in the simulations below). The phase factor is most relevant at places where both E * m (r, φ) and E m (r, φ) are maximized. Therefore, inaccuracies in the center of the disks are in general most relevant.
We parametrize ∆z as a random function where σ is the root-mean-square of the elevation and ξ the transverse correlation length, i.e., the standard deviation characterizing the radius of a typical bump, cf. where ... denotes the average and P 0 is the mixing matrix for the unperturbed system. Thus in the limit of large ξ there is only an overall phase error analogous to a misplacement of the disks. On the other hand, for small ξ k −1 c,m , one finds where the system is dominated by an effective loss, while the phase errors are averaged out. This effective loss can be parameterized within the disks as (analogously to the definition of δ d in eq. (2.7)) with the disk thickness d. The equation holds analogously within the free space gaps, but there it gives around two orders of magnitude smaller δ ∆z (because of = 1 and d ≈ cm there). The effects in both of these limits for ξ can be estimated using 1D calculations as in [22]. For the intermediate range where ξ ≈ k −1 c,m both phase errors and effective loss are relevant. In addition, P will not be well approximated by a diagonal matrix anymore ('mode mixing'), which gives the strongest design constraints, as we will see below. We evaluate representative elements of the mixing matrix for different correlation lengths explicitly in appendix C.
We extend these estimates with explicit numerical results shown in figure 4. We survey the effect on the benchmark power boost factor for both the MADMAX prototype (left column) and full-scale MADMAX (right column). To this end we consider uniformly distributed random tilts γ of the dielectric disks around both x and y axis through their center (first row), non-planar disks (ξ = ø/10, second row) and surface roughness (ξ ≈ λ/4, third row). For each case we take many random samples of a respectively deformed booster and calculate the boost factor times antenna coupling for a Gaussian beam antenna as discussed above.
Calculations  In order to leave the power boost factor β 2 unchanged on the level of 20 % for each individual effect considered alone, we conclude that tilts at the order of γ 100 µm/ø are required, planarity on length scales of ξ ≈ ø/10 should be σ 5 µm and surface roughness is allowed to be σ 20 µm. In an engineering context often the deviation between the minimum and maximum (min-to-max) ∆z min−max is quoted instead of σ. We note that for planarity ∆z min−max ≈ 4σ, and surface roughness ∆z min−max ≈ 6σ.
By defining the tilt around a central axis we have ∆z = 0 and thus suppressed phase errors for the tilts here. This separates the requirement on the tilt from the overall position accuracy of the dielectric disks, which gives more stringent constraints ( 5 µm in this case [52,62]). The strongest design constraints in this section arise for planarity. This is intuitive, since the considered transverse thickness variations ∆z appear on similar length scales ξ as the most relevant modes, maximizing mode mixing effects. The results for surface roughness are consistent with 1D calculations taking losses at the order estimated in eq. (4.6) into account. These constraints remain approximately unchanged when increasing the number of dielectric disks from 20 to 80 disks but keeping the desired boost factor bandwidth the same. This is expected, since a bandwidth of 50 MHz naively corresponds to a resonance with the beam experiencing about 20 GHz/50 MHz ≈ 400 bounces before leaving the booster independently of how many disks are actually installed.
These results are expected to generalize to boost factors at other frequencies with the same relative boost factor bandwidth (ν/∆ν ≈ 400) when written in units of the wavelength λ. This can for example be seen from eq. (4.4), which gives the same P at different frequencies when scaling ∆z accordingly. The above constraints then read: maximum tilts at the order of γ 7 × 10 −3 λ/ø, planarity of σ 4 × 10 −4 λ (on length scales of ξ ≈ ø/10) and maximum surface roughness of σ 1.5 × 10 −3 λ.
These results hold when considering each effect alone. Since the deformations at different length scales are statistically independent, the systematic uncertainty in the boost factor will approximately add in quadrature. Hence, when combining the above constraints, they are expected to tighten by a factor of around √ 3 ≈ 1.7. On the other hand, in an experimental setup one would for example measure the reflectivity in order to constrain the boost factor. Such measurement can be used to realign (tune) the dielectric disks to more optimal positions [50]. Preliminary calculations show that this could approximately soften the planarity constraints by a factor of 2.
Besides the effect on the power boost factor we show the impact on the beam shape of these effects in figure 5. It is evident that the different deformations alter the beam shapes on a similar scale as the size of distortions of the disks in the booster as expected.

Tiled Dielectric Disks
In order to achieve dielectric disks with a diameter of 1 m and low loss, the MADMAX collaboration is also investigating the possibility of gluing together smaller hexagonal patches of LaAlO 3 ( ≈ 24) wavers with a diameter of around 5 cm [24,59]. The gaps between the tiles are filled with glue ( ≈ 5, similar to Stycast 2850FT [67]) and have a thickness of around 0.2 mm, cf. figure 6 (left). In this section we present a first study of the impact of such a tiling on the prototype and full-scale MADMAX benchmark boost factors.
The large difference between dielectric constants on short scales across the glued gaps invalidates the assumption of zero net charge and leads to polarization effects as already seen in [50]. In order to apply the formalism described above, we need to derive a set of eigenmodes of the tiled disks. This can be done semi-analytically with a transfer matrix formalism [56],   Full Scale MADMAX (9 tiles) Ideal Disks Concentrically Tiled Disks but tends to become numerically unstable due to, again, the large relative differences in dielectric constants. This is outside the scope of this work and is left for future studies. However, we can efficiently simulate an azimuthally symmetric geometry with the 2D3D FEM approach introduced in [50]. Therefore, here we consider azimuthally symmetric, concentric tiles as shown in figure 6 (right). The parameter r gap describes the radial distance between two tiles and the gap thickness between two tiles is set to d gap = 0.2 mm. For the prototype we set r gap = 4 cm to approximate the structure shown in figure 6 (left) (corresponding to three gluing gaps for ø = 30 cm, the outermost tile has a width of ≈ 5 cm). For the full-scale MADMAX setup we assume disks with r gap = 6 cm, corresponding roughly to the largest possible diameter of LaAlO 3 crystals with currently available crystal growing techniques [68] (eight gluing gaps for ø = 1 m, width of the outermost tile ≈ 5 cm). Figure 7 shows the result of this calculation for the prototype (left) and full-scale (right) MADMAX benchmark boost factors analogously to figure 2. First, we see that the achievable power boost is only mildly reduced at the level of a few percent compared to the ideal 3D calculation in terms of total emitted power. In addition, the boost factors of the untiled and tiled systems are shifted against each other in frequency. This is consistent with the expectation of additional transverse momentum to the electromagnetic wave obtained from the tiling structure. The shift is much smaller than in the case where each tile would be totally electromagnetically decoupled from each other. In this case the shift according to eq. (3.1) would naively increase by a factor of (ø/r gap ) 2 ≈ 60 (for the prototype) and ≈ 300 (for the full-scale experiment). Next, we consider the emitted beam shape. For the MADMAX prototype the emitted power can still be received with a high efficiency of > 90 % using the Gaussian beam antenna discussed above. However, for the full scale MADMAX we find that the beam shape is significantly altered due to polarization effects caused by the tiling. This is demonstrated in figure 8 where we show the emitted beam shapes of the final scale MADMAX booster at representative frequencies for the full-scale MADMAX boost factor as in figure 7 (right). The electric fields have a non-negligible x-component, although the external magnetic field is polarized in y-direction, B e ∝ê y . At the lower frequency the beam shape is approximately proportional to cos φê φ , at the higher frequency it is dominated by a sin φê r component. At intermediate frequencies it contains both polarizations but at arbitrary phase. Adding them in phase would give a field polarized in y-direction as desired, i.e., contributions from both r and φ polarizations appear shifted with respect to each other in frequency. The x-component always obeys a quadruple structure as we have already seen in [50], cf. also appendix D. Fields ∝ê r are orthogonal and ∝ê φ parallel to the glue gaps and hence need to obey different electromagnetic boundary conditions. Our observations are consistent with r and φ polarized waves therefore having different propagation constants within the booster analogous to propagation in e.g. anisotropic media.
Using the same Gaussian antenna as proposed in the previous sections would therefore reduce the antenna coupling by roughly a factor of 2 or more depending on frequency for the considered full scale MADMAX setup. However, this reduction may be mitigated to some extent by optimization of the antenna shape or disk tiling geometry, as polarization effects may be reduced by using the proper orientation of gaps and shape of tiles. Since resolving the modes for a tiled disk is numerically challenging as described above, studies on alternative tiling designs are not presented here and left for future work.
In addition, reducing the gap size can reduce tiling effects. Also, tiling effects are reduced when decreasing the difference of dielectric constants between the glue and the disks by, for example, using a higher-glue. Alternatively, it is possible to resort to dielectrics with lower dielectric constant , which can also be grown to larger diameters. For example, using sapphire ( ≈ 9) instead of lanthanum aluminate ( ≈ 24) could thus be a realistic alternative to circumvent significant tiling effects, while reducing the power boost factor by acceptable ≈ 30 %. We also note that other communities have experience in building meterscale telescope lenses with high accuracy and without the need of tiling [69]. MADMAX could potentially compensate a reduction in β 2 to some extent by using proportionally more dielectric disks, corresponding to a respective increase in axion-photon conversion volume [22]. Finally, it is noted that MADMAX sensitivity estimates, e.g. in [24], use a conservative system noise temperature of 8 K, which can likely be improved for example by using traveling wave parametric amplifiers [59,70] and would allow for using smaller boost factors.
The above studies only present first estimates for a specific case of tiling. More detailed studies are underway to understand the modes of a tiled booster and the dependencies on the tiling design (such as glue thickness, orientation of tiling gaps, etc.) but also on frequency, on boost factor bandwidth, on the disk diameter and on other parameters. They will provide a clearer picture of the optimal disk design for the full-scale MADMAX booster.

Summary and Conclusion
In this paper we have studied 3D effects in dielectric haloscopes in terms of independently propagating booster eigenmodes. We have derived expected beam shapes for the MADMAX dielectric haloscopes for the first time. The electromagnetic fields inside the booster are not well described by a plane wave, as in previous 1D calculations. However, for finite sized, isotropic and perfectly flat disks the dominant contribution is from the fundamental mode that has a coupling efficiency of 69 % to the axion-induced electric field. This mode can be well received using Gaussian beam quasi-optics [58] matched to a Gaussian beam with a beam waist radius of w 0 ≈ ø/3 focused at the front-most disk of the booster. This can be achieved by using a Gaussian beam horn antenna and one or more focusing mirrors, see e.g. [24,58,59] Table 2: Summary of requirements for MADMAX dielectric haloscopes derived in this paper such as to leave the benchmark boost factor (50 MHz bandwidth at 22 GHz corresponding to m a ≈ 90 µeV) unchanged at the level of 20 % or below.
Moreover, we have derived analytical expressions to quantify the impact on dielectric haloscopes from non-ideal effects such as non-zero axion velocity, magnetic field inhomogeneity and geometrical disk inaccuracies. We have also deduced explicit requirements for the MADMAX dielectric haloscopes for a benchmark boost factor at 22 GHz (m a ≈ 90 µeV) and bandwidth of 50 MHz. Table 2 summarizes these parameters for both the MADMAX prototype and the full-scale MADMAX booster. All values corresponding to the non-ideal booster reduce the boost factor by less than 20 % compared to the ideal 3D case. Realistic values for axion velocities and magnetic field inhomogeneities are mostly unproblematic for MADMAX. However, geometrical inaccuracies of the dielectric disks, such as tilts, non-planarities and surface roughness, cause phase errors, mode mixing and effective losses, and therefore lead to important design constraints. For fixed relative boost factor bandwidth corresponding to 50 MHz at 22 GHz, the results hold approximately independent of the disk number in the considered range between 20 and 80 disks, scale with disk diameter ø as indicated and can be generalized to other frequencies, i.e., axion masses, by appropriate scaling with the wavelength λ. We also have shown that concentric tiling does not reduce the boost factor significantly but shifts it to higher frequencies and can affect the beam shape. Future studies will incorporate polarization effects caused by anisotropic dielectric constants and more realistic tiling designs, such as hexagonally tiled disks.
These results are of significant importance for the experimental design of the MADMAX booster and provide first quantitative design goals for booster manufacturing. In addition, the results may be applicable to other dielectric haloscopes and similar setups such as OR-PHEUS [71,72], DALI [73] and LAMPOST [29].

A.1 Ideal 3D Booster
In order to verify our numerical methods against each other, we have compared the result from Mode Matching with the corresponding results from the Recursive Fourier Propagation and 2D3D FEM methods introduced in [50]. Figure 9 shows this comparison for the MADMAX prototype benchmark boost factor discussed in this paper. The 2D3D FEM method solves the full vectorized wave eq. (2.1), while the other methods assume a scalar diffraction theory and neglect free charges in the following by setting ∇ · E = 0 here. The Mode Matching method in addition neglects higher modes, here m > 5, > 2. The lower panel shows the relative difference between the results from Recursive Fourier Propagation and 2D3D FEM against the result from the Mode Matching method, while in the upper panel the boost factor obtained with Mode Matching is shown for orientation. The systematic differences are likely due to the simplifying assumptions of the Recursive Fourier and Mode Matching methods. Most prominently, the boost factors obtained by Recursive Fourier Propagation and FEM are typically higher than results from Mode Matching. This is expected, since the Mode Matching method neglects higher modes which may carry additional power. The differences are largest in the regions where the boost factor itself is small. The boost factors are consistent up to percent level within the boost factor bandwidth. This does not significantly affect sensitivity and therefore is sufficient for this study.
It should be noted that the largest deviations indeed are outside the 50 MHz range of the boost factor, where higher modes contribute. Analogous results have been obtained for the MADMAX prototype within its designated frequency range at ν = (18,20,22,24) GHz and at 22 GHz for different boost factor bandwidths of (5, 10, 20, 50, 100, 250) MHz. In addition, the comparison has also been performed for the full-scale MADMAX setup for the 50 MHz benchmark boost factor at 22 GHz shown in figure 2 (right), leading to analogous results.
The agreement shows for the ideal 3D dielectric haloscopes consisting of multiple finitesized disks tuned to a boost factor over a bandwidth 10 −3 ν the simplifying physics assumptions of the Mode Matching and Recursive Fourier Propagation methods are valid, i.e., a relatively low number of modes (in our case 4) is sufficient to approximate the fields inside the system and a scalar diffraction theory neglecting effects from free charges is sufficient.

A.2 Non-Ideal Booster
The calculations presented in section 4 are not feasible with the 2D3D FEM method, since the azimuthal symmetry is broken for these non-ideal boosters. Therefore, explicit numerical confirmation of the scalar diffraction theory for the non-ideal booster remains for future work. However, the scalar theory holds in the limit where k ⊥ E, which is still a good approximation for the MADMAX setups discussed here. In addition, we have compared results from Mode Matching with results from Recursive Fourier Propagation for the 20 disk benchmark boost factor as in figure 2 (left). Figure 10 shows the same result for an exaggerated pessimistic planarity of σ = 10 µm at a scale of ξ = 35 mm. The percent level differences are irrelevant for sensitivity estimates. We also show a comparison for the planarity calculation between the beam shapes obtained at the frequency with maximum boost in figure 11. The observed differences are at smaller scales than the considered modes, i.e., mainly arise due to the fact that the Mode Matching method is neglecting higher modes. The differences can be reduced when taking into account more modes. Analogous results are obtained for velocity effects, magnetic field inhomogeneities and tilts.

C Explicit Mode Mixing Matrix Calculations for Transverse Disk Thickness Variations
We have verified the initial estimates in section 4.3 by calculating the mixing matrix P numerically for many realization of thickness variations ∆z(r, φ) and observing how it changes with ξ. Figure 12 shows the result of such a calculation for σ = 10 µm. Other σ give analogous results. Here we used ≈

D Allowed Fields for Azimuthal Symmetry and Linearly Polarized Source Term
In the case of an azimuthally symmetric geometry and a linearly polarized external B-field in the y-direction, we can use the 2D3D approach introduced in [50]. The total solution is obtained as a superposition where the fieldsẼ + andẼ − are calculated numerically as described in [50]. It is important to notice that the ± solutions obey the relations which have not been used in [50]. Using eq. (D.1)-(D.4) theê r ,ê φ andê z contributions to the total E are From eq. (D.5)-(D.7) we see that the r-polarized electric fields always have sin φ and the φ-polarized ones cos φ dependence. We also note that the x-component of the fields obeys always a quadrupole structure. The y-component of the r (φ) component obeys a vertical (horizontal) dipole structure.