1 Introduction↩︎

The nature of DM, constituting most of the mass in our Universe, remains an elusive issue [1]. Through more than three decades, the most promising candidates have been considered to be Weakly Interacting Massive Particles (WIMPs), with masses in the GeV-to-TeV range. Yet, after decades of searching, no conclusive detection of WIMPs has been made [2]–[4]. As a consequence, the search for DM has diversified substantially in the last decade [5], with lighter DM particles gaining much attention, both theoretically and experimentally [6].

The direct detection of sub-GeV DM particles poses major experimental challenges. Standard searches for DM-nucleus scattering become ineffective at low mass, as the low DM velocity in the Milky Way Halo (\(\mathcal{O}(10^{-3})\,c)\)) [7] and kinematic mismatch between the DM and nuclear masses limit the amount of energy which can be deposited. In recent years, several detection channels for light and ultralight DM have been theoretically proposed, such as the excitation of phonons, polaritons or magnons [8]–[10] or the absorption of DM in semiconductors [11], [12] and superconductors [13], [14]. In general, collective excitations must be taken into account when the de Broglie wavelength becomes larger that the particle spacing; when dealing with solid targets for DM detection, this corresponds roughly to particle masses below 1 MeV, becoming progressively relevant as the particle mass decreases. Furthermore, collective excitations can additionally imply a daily modulation in the signal if anisotropic targets are used [15]. This could provide further evidence for DM, allowing us to more easily distinguish a genuine DM signal from unmodulated background noise.

Among newly proposed strategies for the direct detection of low mass DM, of special interest is the excitation of optical phonons in polar materials [8], because of their significant coupling to light dark photon mediators. Moreover, optical phonons in polar materials display gaps of the order of 10-30 meV, favourable for the scattering of DM particles with masses below 1 MeV. The detection of such low energy excitations requires extremely high sensitivity, only achievable by cryogenic detectors. Indeed, cryogenic thermal sensors, operating at temperatures well below 1K and using a phonon signal to determine the total energy released in an absorber/target, have become essential to build up very sensitive rare event detectors for fundamental physics experiments and astrophysics [16]. Transition Edge Sensors (TES) [17] are widely used in or proposed for telescopes both on Earth and in space, to detect photons over a broad range of the electromagnetic spectrum. TES have also been used for years - or are being proposed - for particle detectors, such as those for dark matter (DM) searches (CRESST [18], [19] , SuperCDMS [20], [21], COSINUS [22]) or neutrino physics (HOLMES [23], PTOLEMY [24]). Other types of cryogenic detectors are also used or in development, in a variety of initiatives: NTDs (EDELWEISS [25]); kinetic inductance detectors (KIDs) (BULLKID [26], CADEx [27], CALDER [28]); superconducting single photon nanowires for light DM [29]; or magnetic microcalorimeters for neutrinos (ECHo [30]). Many cryogenic particle detectors are based on the so-called quasi-particle trap-assisted TES, QET [31] (Figure 1), in which a superconducting film, placed between the target and the TES, is used as a phonon collector. Its role is to decouple the phonon collection efficiency (related to the target surface area covering) from the TES sensitivity (proportional to TES volume), thus allowing independent optimization of both.

The detection of single phonons resulting from sub-MeV DM scattering will require fine tuning of detector performances (see e.g. Fink and coworkers [32]), involving understanding and optimization of all the steps of the detection process. While the reach and daily modulation of the signal are mostly determined by the harmonic properties of the target material, the ultimate performances of the detector, in terms of efficiency and phonon energy sensitivity, are limited by a number of factors. These include the TES sensitivity and time response (\(t_r \gtrsim \mu\mathrm{s}\)), as well as the processes undergone by the phonons initially generated in the target as a result of the DM interaction (sometimes referred to as athermal phonons, to stress that they are strongly out-of-equilibrium quasi-particles that do not correspond to the equilibrium thermodynamic conditions of the target material). This means that knowing the phonon properties of the target (spectra, transport, role of defects, interface transmission rate) is crucial for designing and building a detector capable of detecting low mass dark matter through single phonon excitations.

A benchmarking of polar materials as possible targets has been published recently by Griffin and coworkers [33], taking into account different detection channels. There, the expected detector reach has been evaluated considering several types of DM candidates and using density-functional theory (DFT) calculations of phonon properties. This initial work helps in identifying the material parameters that determine the detector reach, and to make a pre-selection of target candidates. Calculations at the harmonic level (i.e.phonon band structure) are indeed very useful as a pre-screening tool, but other relevant features can only be addressed by means of anharmonic calculations, where the cross-section of the relevant scattering processes are computed [34]. In particular, the mechanisms for phonon downconversion, the conditions for the subsequent ballistic transport of long wavelength phonons and the role of isotope purity cannot be tackled at the harmonic level and have thus far been neglected. Only a recent paper [35] explores the phonon transport in Si with GEANT [36], a software tool used in high energy physics to study particle-matter interactions. Although they provide some useful indications, these simulations are based on a number of oversimplifying assumptions that call into question the predictive power of these results. Similarly, phonon reflection at the interface with the superconducting collector is often neglected or treated within simplified phenomenological models.

In this paper we study theoretically, within an entirely ab initio scheme, the processes undergone by phonons created by the interaction of a light DM particle in a polar semiconductor (target), until they reach the phonon collector. These results are essential to assess the feasibility of an innovative light DM detector and to evaluate (and then minimize) its non-ideal performances, i.e.the loss of phonon signal reaching the sensors. We focus on Al\(_2\)O\(_3\)which, as we discuss below, checks all the boxes of the ideal polar target for DM detection. The phonon collector is taken to be Al, which, besides its suitable superconducting properties, presents high-quality interfaces with Al\(_2\)O\(_3\). However, we emphasise that we focus exclusively on the phonon dynamics within the target, and that the properties of the collector only enter in the determination of the phonon transmission probability across the Al\(_2\)O\(_3\)/Al interface.

The structure of the paper is the following. In Section 2 we give an overview of DM phonon scattering mechanisms and the detection scheme. In Section 3 we present and discuss our results; specifically, (i) we provide an exhaustive ab initio characterization of the phononic properties of the target/collector system chosen, i.e. Al\(_2\)O\(_3\)/Al; (ii) by time integration of the collision operator we discuss the validity of a number of approximations that can be adopted to describe the phonon downconversion and subsequent propagation. Finally, we present results from a newly developed energy-deviational Monte Carlo scheme, assessing the performances of the detector and discussing the feasibility of the detection strategy. In Section 4 we present our conclusions. The appendices provide a more insightful discussion regarding the employed methodologies.

2 Phonon-based direct detection of light dark matter↩︎

The production of phonons offers an attractive channel for searching for light Dark Matter particles, especially for DM particles lighter than \(\sim 1 \,\mathrm{MeV}\). The de Broglie wavelength of such light particles is much larger than the inter-atomic spacing in crystal targets, meaning that they can effectively couple with phonons in the material [37]–[39]. While sub-MeV particles typically do not have sufficient kinetic energy to excite detectable nuclear recoils or electronic ionization, the production of phonons (with energies \(\omega \gtrsim \mathrm{meV}\)) is instead kinematically accessible.

Depending on the model of DM, phonons may be produced by scattering or by the absorption of the DM particle directly by the target material. The rate of production of phonons with energy \(\omega\) and momentum \(\mathbf{q}\) can be obtained by convolving the DM-phonon interaction rate \(\mathrm{d}\Gamma/\mathrm{d}\omega\mathrm{d}\mathbf{k}\) with the velocity distribution of DM in the Milky Way halo \(f(\mathbf{v})\). Writing this rate per unit time and per unit target mass, we obtain [15]: \[\label{eq:DMRate} \frac{\mathrm{d}R}{\mathrm{d}\omega\,\mathrm{d}\mathbf{q}} = \frac{1}{\rho_\mathrm{T}}\frac{\rho_\mathrm{DM}}{m_\chi} \int \mathrm{d}^3\mathbf{v} f(\mathbf{v}) \frac{\mathrm{d}\Gamma (\mathbf{v})}{\mathrm{d}\omega\,\mathrm{d}\mathbf{k}}\,.\tag{1}\] Here, \(\rho_\mathrm{T}\) is the mass density of the target material. The local density of DM can be estimated from mass modeling of the Milky Way or from the dynamics of local stars [40], [41], giving \(\rho_\mathrm{DM} \approx 0.3 - 0.6 \,\mathrm{GeV}\,\mathrm{cm}^{-3}\). The DM velocity distribution in the Earth frame is typically assumed to take the form of a Maxwell-Boltzmann distribution (with typical velocity \(v_0 \approx 220 \,\mathrm{km/s}\)) and a cut-off at the Galactic escape velocity (\(\sim 775\,\mathrm{km/s}\) in the Earth frame) [42]. A number of public tools have been made available for the numerical calculation of the rates described by Eq. 1 , including PhonoDark [43] and DarkELF [44].

As a concrete example, we will consider the case of a DM particle interacting via an ultra-light dark photon mediator. This effectively gives the DM particle a very small electromagnetic charge. Calculation of the DM-phonon scattering rate then proceeds in analogy with that of electron-phonon scattering, rescaled by the small effective charge of the DM [8]. This scattering rate can also be related to the energy loss function (ELF) of the target material, which is well-studied both theoretically and experimentally [45], [46]. One advantage of DM particles interacting via a dark photon is that they preferentially couple to optical phonons, which is particularly promising. For a fixed DM mass and velocity, optical phonons typically allow for a larger energy deposition compared to acoustic phonons, which comes purely from kinematic considerations and the fact that the optical phonons are gapped [38]. We therefore focus here on the production of a single optical phonon from DM interactions. The scattering of DM to produce two optical phonons is expected to be subdominant, and even kinematically inaccessible for sufficiently small DM masses [47]–[49].

In order to detect the athermal phonon produced in the DM interaction, its energy must be transported to the surface of the target crystal and transferred to the phonon collector. The overall process can be broken down in different stages (see the sketch in Figure 1). At first the optical phonon generated by the DM particle must downconvert to a lower energy, long wavelength acoustic phonon. This is necessary for two reasons: (i) high-energy optical phonons have typically very small group velocities and do not propagate efficiently; (ii) if the phonon energy does not decrease, the elastic mismatch with the collector prevents its transmission. After downconversion, the resulting long wavelength phonon travels quasi-ballistically through the target (something made possible by the very low operation temperatures) until it reaches the collector. There, the phonon can be transmitted or reflected, with a probability that depends on the thermal boundary resistance between the target and the collector. If it is reflected, it will cross the whole length of the target twice (once backward, once forward after reflection from the backsurface) and attempt again to be transmitted. Typically, these multiple reflections occur when the initial downconversion process did not lower sufficiently the phonon energy, and the additional travelling through the target offers more opportunities for further downconversion events. The transmission probability at the interface is essentially ruled by mismatch of the elastic properties of the target and the collector, as at such low temperatures anharmonic effects can be likely neglected.

In the following, we study these phonon processes in more detail, beginning with a description of the phonon properties of Al\(_2\)O\(_3\)and Al, which we consider as promising materials for the target and phonon collector, respectively.

3 Results and Discussion↩︎

3.1 Materials: phonons and phonon-related properties for DM detection↩︎

Phonon dispersions. The starting point to approach any phonon-based detection scheme is the phonon dispersion of the materials involved. Despite the increasing popularity of large-scale inorganic material databases, only a few of them include information about phonon properties (e.g. Materials Project [50], almaBTE [51], phonondb@kyoto-u [52]). Here we compute the phonon dispersion of Al\(_2\)O\(_3\)and Al by finite differences in \(4 \times 4 \times 4\) and \(8 \times 8 \times 8\) supercells, respectively.

We use the code [54] to generate all the inequivalent atomic displacements (in such supercells) that are required to compute the second order interatomic force constants (IFC) and then the dynamical matrix. The results for Al\(_2\)O\(_3\)and Al are shown in Figure 2. A few things are worth noting: (i) the bandwidth of the phonon dispersion of Al\(_2\)O\(_3\)covers the energy range that is expected to be important for sub-MeV DM detection, i.e.\(\omega \lesssim 100 \,\mathrm{meV}\); (ii) the target is anisotropic (e.g. the phonon bands along \(\Gamma-Z\), \(\Gamma-A\), and \(\Gamma-D\) are different) and is thus suited to produce a daily modulation of the signal due to DM particles, which would provide further evidence of their detection; (iii) the spectra of Al\(_2\)O\(_3\)and Al feature a considerable energy mismatch (phonons in Al cover barely one third of the bandwidth of Al\(_2\)O\(_3\)). This mismatch will play an important role when it comes to compute the transmission probability of phonons when they hit the target/collector interface (see below).

It is important to highlight that having phonon states in the right energy range for the detection of light and ultralight DM is only a minimal requirement. In reality, one should convolve the DM-phonon scattering rate distribution of the DM particles incident on the target, as described in Eq. 1 . As a concrete example, we show in Figure 3, the expected rate of production \(\mathrm{d}R/\mathrm{d}\omega\mathrm{d}k\) of phonons with energy \(\omega\) and momentum \(k\), due to DM-phonon scattering. We perform the calculation using the code [46], assuming a fermionic DM particle of mass \(m_\chi = 100 \,\mathrm{keV}\) interacting via an ultralight dark photon mediator. For simplicity, we assume phonon production only along the ordinary direction (i.e.the \(\langle 0001 \rangle\) direction) of Al\(_2\)O\(_3\).

At a given phonon energy, the rate of phonon production is non-zero only for a finite range of momenta. This range is set by kinematics; requiring conservation of energy and momentum in the scattering process leads to the constraint \(\omega \leq \lVert \mathbf{q}\rVert v_\chi - \lVert \mathbf{q}\rVert^2/2m_\chi\), where \(v_\chi\) is the incoming DM velocity. The maximum value of \(v_\chi\) is set by the sum of the Earth’s velocity and the Milky Way escape velocity (\(v_\mathrm{max} \approx 775\,\mathrm{km/s}\)) [7], [42], and this in turn sets the maximum momentum for phonon states which can be produced in the scattering event. Overplotted as red dots in Figure 3 are the available phonon states in the \(\mathbf{q}\)-mesh we use in this work, demonstrating that some of these are accessible.

Thermal conductivity. The detection process begins with a DM particle scattering in the target to produce an optical phonon. For the detection to be efficient, we need this phonon first to downconvert to longer wavelength, quasi-ballistic phonons, which, in turn, should travel long distances with minimal scattering, so that a large signal reaches the target/collector interface. These two requirements are to some extent contradictory: downconversion relies on a certain level of anharmonicity, but anharmonicity is also the main cause for phonon scattering. To have a first rough idea, we computed the thermal conductivity of Al\(_2\)O\(_3\)by solving the Peierls-Boltzmann Transport Equation (PBTE) [34], [55], [56], using as inputs the harmonic and anharmonic interatomic force constants (IFCs) obtained from DFT calculations. The third-order anharmonic IFCs were calculated in the same supercells used for the harmonic ones, considering interactions up to 7\(^\text{th}\) neighbors. Scattering from isotopes is accounted for within the model of Tamura [57], considering the natural isotopic population of O [Al has only one stable isotope, \(^{27}\)Al, while O has three: \(^{16}\)O (99.8%), \(^{17}\)O (0.038%), and \(^{18}\)O (0.205%)]. We solve the PBTE on a \(17 \times 17 \times 17\) q-point grid with the almaBTE code [58] (see Methods for full details of the calculations).

We obtain a thermal conductivity at 300 K of \(\kappa_{xx}=\kappa_{yy}=44.0\) and \(\kappa_{zz}=48.2\) W m\(^{-1}\) K\(^{-1}\), which increases up to 500-550 W m\(^{-1}\) K\(^{-1}\) when the temperature decreases down to 100 K, as phonon-phonon scattering processes become much weaker (see Figure 4). These results agree well with those previously obtained by Dongre and coworkers [59]. Clearly, at the temperatures relevant for the present detection scheme, \(T \approx 10\) mK, the thermal conductivity will be much higher. However, a reliable estimate of \({\boldsymbol{\kappa}}\) at such low temperatures is beyond current computational capabilities, because of the extraordinarily fine q-point mesh necessary to properly capture very long wavelength phonons, which become increasingly important as the temperature decreases.

These results simply provide a general indication that Al\(_2\)O\(_3\)has a medium room temperature thermal conductivity for a single-crystalline solid. What matters for the detection scheme under scrutiny is how a specific optical phonon downconverts and how efficiently the (ideally quasi-ballistic) phonons resulting from the downconversion process propagate. These questions will be specifically addressed below.

The insets of Figure 4 report two results that are important for the forthcoming discussion. On the left hand-side we plot the ratio between the values coming from a full solution of the PBTE and those obtained within the simpler Relaxation Time Approximation (RTA) [60] (see Sec. 3.2). Corrections to the RTA made by the full solution of the PBTE are essentially negligible at room temperature, but they increase up to \(\approx 6\)% at 100 K, a temperature which is still four orders of magnitude larger than the operating temperature of a phonon-based DM detector. The inset on the right displays the thermal conductivity anisotropy, \(\kappa_{zz}/\kappa_{xx}\), which, as discussed above, is important when it comes to the daily modulation of the detected signal.

Thermal boundary resistance. When the heat flux hits the target/collector interface, part of it will be transmitted with a certain probability, \(\alpha\), and part will be reflected. The fact that not all phonons are transmitted can be understood in terms of the mismatch of the elastic properties of the materials on the two sides of the interface. This is easy to understand looking at the phonon dispersions of Al\(_2\)O\(_3\)and Al displayed in Figure 2: it is clear that phonons of the target (Al\(_2\)O\(_3\)) with \(f > 10\) THz cannot be transmitted to the collector (Al), because there are no available phonon states at that energy. The situation is further complicated when considering that (i) even when a finite vibrational density of states exists on the two sides of the interface, the transmission probability depends on the degree of coupling of the states involved and is in general lower than unity, and (ii) the effect of interface anharmonic scattering can open additional transmission channels [61]. Indeed the development of a microscopic theory of thermal boundary resistance (TBR) is one of the most active areas in thermal science [62], [63].

In the following, we will make use of one of the most common phenomenological models for the calculation of the TBR, the so-called diffuse mismatch model (DMM)  [62] in a full-band formalism [64], [65]. Within DMM, all phonons undergo fully diffusive scattering at the interface¹ and they lose memory of their origin after being scattered. In Figure 5 we plot the transmission probability across an Al\(_2\)O\(_3\)/Al interface for each phonon on the Al\(_2\)O\(_3\)side. As anticipated, the transmission of phonons with \(f > 10\) THz is suppressed, simply because no states at that energy are available on the Al side. Transmission at lower frequencies is larger, though a finite reflection probability exists in most cases.

3.2 Time integration of the collision operator↩︎

The Peierls-Boltzmann Transport Equation. As the common choice for target materials are high-quality single-crystal polar semiconductors, the phonon dynamics can be in principle decoupled from the electronic one, and thus can be described through the deviational Peierls-Boltzmann Transport Equation (PBTE) [55], [56]: \[\label{PBTE} \frac{\partial n^d_i}{\partial t} + \mathbf{v}_i \cdot {\nabla_r} n^d_i + \frac{\partial n^0_i}{\partial T} \mathbf{v}_i \cdot {\nabla_r} T_{\mathrm{ref}}= \left.\frac{\partial n_i}{\partial t} \right|_\mathrm{collision}\,.\tag{2}\] This approach is deviational in the sense that rather than describing the evolution of the total phonon population in the target, \(n\), which normally results in poor signal-to-noise ratios, it deals with the deviation of the phonon population from equilibrium, \(n_d = n - n_0\). The equilibrium distribution \(n^0\) is well described by a Bose-Einstein distribution at the target temperature \(T_{\mathrm{ref}}\), while the deviation \(n^d\) is the DM-generated (athermal) phonon contribution. Here, \(\mathbf{v}_i\) is the group velocity of the \(i\)-th mode; and \(\left.\frac{\partial n_i}{\partial t} \right|_\mathrm{collision}\) is the collision operator, modeling the change in the phonon distribution due to all the possible phonon scattering mechanisms, e.g.phonon-phonon processes, isotopic and boundary scattering.

Due to its complexity, the full nonlinear phonon collision operator presents significant challenges when solving the PBTE, making its practical use limited in most situations. A common approach consists in reducing it to a more tractable form by linearizing it, i.e.by considering only small deviations from equilibrium, \(n^0 \gg n^d\). The collision operator reduces then to: \[\label{eq:collisionAnd} \frac{\partial n_i^d}{\partial t}\Big|_{\mathrm{collision}} = A_{ij} n_j^d\tag{3}\] where \(A_{ij}\) is the linearized collision operator.

Sometimes, the PBTE is further simplified using the so-called relaxation time approximation (RTA), where it is assumed that each phonon mode relaxes to equilibrium independently with a characteristic time, \(\tau_{i}\). This translates to supposing that only the involved mode, \(i\), is slightly out of equilibrium, while the other modes remain in equilibrium (i.e. \(n_j = n_j^0 \; \forall\; j\neq i\)). In this way, the RTA collision operator, \(A_{ij}^\mathrm{RTA}\), can be expressed in terms of the full linearized collision operator as: \[A_{ij}^\mathrm{RTA} = A_{ij}\delta_{ij} = -\frac{n^d_i}{\tau_i}\delta_{ij},\] where the last is the most common expression for the RTA.

Despite its simplicity, which has earned it a great renown in the field of phonon and lattice thermal calculations (particularly when fully ab initio calculations of the thermal conductivity first became possible [66]) the RTA has several limitations. It is easy to see that, under this approximation, the collision operator does not conserve heat flux to any degree, deeming all the three-phonon processes resistive. Such a flaw causes the RTA to fail for cases in which non-resistive processes, i.e.those that mostly redistribute the phonon population, play a non-negligible role. This is the case of e.g.diamond and 2D materials, for which the RTA has been observed to provide quite a poor description of thermal properties, even at room temperature. However, at sufficiently low temperatures this situation becomes common to all materials and predictions based on the RTA are in general unreliable [66]–[68].

The target semiconductor material used here, Al\(_2\)O\(_3\), is a good example: the error induced by the RTA, when compared to a full solution of the PBTE, is negligible at room temperature, but it already increases to \(\approx 6\)% at \(T=100\) K (see the discussion above). Although we provide a careful validation of the RTA below, these considerations by themselves already suggest that the RTA is most likely bound to yield an unreliable description of a phonon-based detection scheme that is expected to operate at cryogenic temperatures.

There is another limitation of a pure RTA operator which is specific of this kind of setup that is important to discuss. The central assumption of the RTA is that each nonequilibrium (athermal) phonon mode returns to equilibrium independently from the rest. It follows from this approximation that redistribution among the different modes is not permitted and thus no flux of phonons will be generated if the initial state is fully located at a standing mode (i.e. \(\mathbf{v}_\lambda = \mathbf{0}\)). A more realistic approach in our case would be to allow a relaxation of the initial mode \(K\) to other modes, which then thermalize (i.e.allowing only initial state out-coupling). This generalization of the RTA, which we dub pseudo-RTA, serves in principle our purposes and can be expressed in terms of the linearized collision operator as: \[\label{eq:pseudoRTA} A_{ij}^{\mathrm{pseudo-RTA}} = A_{ij} [i = j ~\mathrm{or}~j = K].\tag{4}\]

The linearization of the PBTE is a much more widespread approximation than the RTA and virtually all the codes that calculate the thermal conductivity from ab initio computed IFCs rely on it, even when beyond-RTA solutions are pursued [54], [58], [69]. The validity of this approximation, which is not normally the source of major problems, should nonetheless be addressed in detail in the case we are studying: it is true that we assume that only one phonon mode is excited from the detection of DM, but at the same time is driven to a strongly out-of-equilibrium condition. It is difficult to say which of these two factors dominates and thus the reliability of a linearized solution, as well as of the pseudo-RTA, are carefully assessed below.
Different treatments of the collision operator. As a first approach to the problem of phonon-based DM detection, we present results of the time integration of the collision operator in our Al\(_2\)O\(_3\)/Al system. These calculations serve to discuss the validity of the two approximations discussed above, i.e. the linearization and the RTA, which have been previously used in the literature [70].

Our starting point is a simplistic expression of the heat flux transmitted from the target material to the collector [62]: \[\begin{align} \label{eq:GrossFlux} J _{\mathrm{gross}}(t) & = & \frac{1}{N_{\mathrm{q}}V_{\mathrm{uc}}}\sum_i \alpha_i J_{\mathrm{proj},i} (t) \nonumber \\ &=& \frac{1}{N_{\mathrm{q}}V_{\mathrm{uc}}}\sum_i \alpha_i \hbar\omega_i n^d_i(t) \mathbf{v}_i\cdot \mathbf{u} ~ [\mathbf{v}_i\cdot \mathbf{u} > 0], \end{align}\tag{5}\] where \(N_{\mathrm{q}}\) is the number of q-points, \(V_{\mathrm{uc}}\) is the volume of the unit cell, \(\hbar\) is the reduced Planck constant, \(\omega_i\) is the frequency of the \(i\)-th phonon mode, \(\alpha_i\) is the transmission coefficient of the given mode to the collector material, \(\mathbf{v}_i\) is the group velocity, and \(\mathbf{u}\) is a normalized vector pointing out of the target into the collector; \(n^d(t)\) is the deviational phonon distribution, i.e. the difference between the nonequilibrium and the equilibrium distribution function, \(n - n^0\).

In Figure 6 we plot the time evolution of the projected heat flux, comparing the different approaches used to treat the collision operator. As it can be seen, relying on the pseudo-RTA results in a severe underestimation of the heat flux. This is not surprising because, as discussed above, the RTA considers all phonon collisions to be resistive, while at the very low temperatures here considered the majority of phonon-phonon process simply redistribute momentum, roughly conserving the heat flux. The two beyond-RTA solutions, on the other hand, describe correctly the phonon scattering dynamics and yield larger fluxes, thus increasing the probability that the signal originated from the DM interaction reaches the collector and can be successfully detected.

Remarkably, we find virtually no difference between the solution obtained from the full and the linearized collision operator (lower inset of Figure 6). This is a very important conclusion, not only because the linearized solution is computationally less expensive, but also because this opens up the possibility to use Monte Carlo algorithms (which are considered the gold standard in solving the BTE, but are only applicable when linearization holds). Notice, incidentally, that the results in Figure 6 were obtained assuming a volume of 1000 nm\(^3\), which is a very small value compared to any realistic target sample. Therefore, these conclusions are particularly robust.

Notice that in the beyond-RTA solutions we obtain a build up of the heat flux, i.e.\(J\) increases with time, due to the subsequent energy downconversion of phonons, which leads to the increase in the population of lower frequency modes with higher group velocities. Such a build up is absent in the pseudo-RTA, where phonons are allowed to downconvert only once.

As linearization seems to be a valid approximation, at least in the limit of a single-phonon perturbation signal, this is the framework assumed from here on. Likewise, any RTA-based approach must be discarded.

3.3 Energy-deviational Monte Carlo simulations↩︎

The results presented in Section 3.2 are important because they already capture important physical effects that are relevant for phonon-based detection of DM. Above all, however, they permitted us to assess different approaches to the collision operator, leading us to conclude that it can be safely linearized, while any RTA-derived description at such low temperatures is necessarily inaccurate.

However, the main limitation of the theoretical framework adopted thus far is that the spatial dimension is completely neglected. Therefore, not even the simplest device geometries can be described. In this scheme, a DM particle is detected and the time evolution of the phonons generated in the target can be studied, but the value of the heat flux is the same at any point in space (or, equivalently, the propagation speed is infinite). In this way, not only are nuances of phonon propagation, e.g.multiple reflections, completely unaccounted for, but even a basic question, such as what the optimal target thickness is, remains unanswered.

To circumvent this limitation, we developed a 3D deviational-energy Monte Carlo method that relies on a description of the linearized collision operator beyond the RTA (see Methods), whose accuracy has been validated in the previous section. In this framework, and after linearizing the collision operator, Eq. 2 becomes: \[\frac{\partial n^d_i}{\partial t} + \mathbf{v}_i \nabla n^d_i + \frac{\partial n^0_i}{\partial T} \mathbf{v}_i \cdot \nabla T_\mathrm{ref} = \sum_j{A_{ij}n^d_j}\, . \label{eq:mc-dev}\tag{6}\] Also, we make the transformations \(f_i = \hbar \omega_i n_i\) and \(B_{ij} = \frac{\omega_i}{\omega_j} A_{ij}\) and obtain the final deviational-energy formalism: \[\frac{\partial f^d_i}{\partial t} + \mathbf{v}_i \nabla f^d_i + \frac{\partial f^0_i}{\partial T} \mathbf{v}_i \cdot \nabla T_\mathrm{ref} = \sum_j{B_{ij}f^d_j} \, , \label{eq:mc-dev-ene}\tag{7}\] which naturally guarantees energy conservation, without the need for additional algorithms that might bias the solution.

The device setup that we study is illustrated in Figure 7. The first question that we address is to what extent the energy and momentum of the optical phonon created as a result of the DM scattering event matters. To this end we selected phonons of Al\(_2\)O\(_3\)that, according to the maps of Figure 3, can be excited by DM and we studied their evolution. Specifically, we have selected the following highly interacting phonon modes: high-energy-low-momentum (\(E=\SI{0.11}{\electronvolt},~\lVert\mathbf{q}\rVert=\SI{171}{\electronvolt}\)), high-energy-high-momentum (\(E=\SI{0.11}{\electronvolt},~\lVert\mathbf{q}\rVert=\SI{342}{\electronvolt}\)), low-energy-low-momentum (\(E=\SI{0.06}{\electronvolt},~\lVert\mathbf{q}\rVert=\SI{171}{\electronvolt}\)), and high-energy-high-momentum (\(E=\SI{0.06}{\electronvolt},~\lVert\mathbf{q}\rVert=\SI{342}{\electronvolt}\)) modes, all of them parallel to the \(\langle0001\rangle\) crystalline direction.

The results are plotted in Figure 8 for two different device sizes, \(L\), and indicate that higher energies provide higher signal, i.e.larger heat fluxes and thus power. Conversely, phonon momentum seems not to be important at all. Indeed, we considered the time evolution of two phonon states of essentially the same energy, \(\omega \sim 0.11\) eV, but different q-vectors, finding that the differences in the generated power are negligible (\(L = 10\) nm, black and red curves in the left panel of Figure 8) or small (\(L = 100\) nm, green and blue curves in the right panel). Notice that we focus our analysis on the tail of the signals, considering the typical response time of TES. Signal peaks, which are not significative for our purposes due to the aforementioned reasons, can be appreciated in the log-log plots in the insets.

Next, now that the spatial dimension can be directly included in our simulations, we study the effect of the location where the DM generates an optical phonon in the target, \(r_0\) (or, equivalently, the distance from the collector, \(d = L - r_0\)). In Figure 9, where we plot the power as a function of time right after the target/collector interface, we found that the onset of the signal trivially depends on the distance that the phonons must travel before they reach the interface, as can be clearly seen in the inset. On the other hand, these differences tend to blur out at long simulation times, when phonons are transmitted after multiple reflections within the target (i.e., between the target/collector interface and the target backsurface) and then small differences in the initial conditions are less important. Nevertheless, the difference tends to decrease at long simulation times, thicker devices result in larger signal (e.g.the power in the case of \(L=1000\) nm is always larger than the one for \(L=10\) or 100 nm). This last can be seen as the direct consequence of the dominance of boundary scattering in thinner devices, which as harmonic, do not allow for decaying to transmissible states, thus reducing the overall collected signal.

We recall that, due to the slow response time of TES, this is the important region for our analysis. Notice that in these calculations we assumed Al to be oriented along the \(\langle 111 \rangle\) axis (see also the sketch in Figure 7); this matters when it comes to calculate the TBR. To make sure that no dramatic changes arise for different interfaces, we have repeated the calculations considering the \(\langle 110 \rangle\) crystallographic orientation for the Al collector, finding that changes are negligible (see Supplementary Figure 1).

Finally, we present what is perhaps the most striking result of our calculations. The common wisdom is that the mass disorder introduced by a given isotope distribution is a source of scattering and should result in lower heat fluxes and thus power. Indeed, one of the reasons for selecting Al\(_2\)O\(_3\)is its small isotope dispersion (Al has only one stable isotope, while 99.8% O atoms are \(^{16}\)O) and a valid question is if one should not rather pick a naturally isotopically-pure semiconductor, e.g.AlP, as a target material or think about isotope enrichment. For instance, Harrelson and coworkers [70] found that isotopic scattering in GaAs and Si-based targets is the dominant scattering process, with associated lifetimes that are significantly smaller than anharmonic ones. They concluded that this effect is mitigated only with isotopic enrichment greater than 99.9%.

At variance with these considerations, however, we found that the main role of isotopes is providing, via scattering, transport channels of higher velocities. This is shown in Figure 10, where we compare the time evolution of the heat power in natural and isotopically purified Al\(_2\)O\(_3\) (i.e.all O atoms are \(^{16}\)O). As it can be seen, the behavior is qualitatively similar, regardless of the values of the energy and momentum of the initially excited phonon: switching off isotope scattering decreases, rather than increases, the heat flux. A reliable estimate of this decrease at large values of \(t\) is difficult to make, because the signal of the isotopically-pure case is so low that it is comparable to MC statistical noise. However, we conservatively estimate such a reduction to be of around two orders of magnitude for \(t \approx 100\) ns. It is difficult to generalize this result to other materials, as it ultimately depends on the isotope population, isotope scattering rates and phonon velocity distribution. Nevertheless, at least in the case of Al\(_2\)O\(_3\), rather than a hurdle, isotopes can even become an ally when it comes to maximizing the signal that reaches the superconducting film.

Daily modulation and phonon dynamics. As Earth rotates the alignment of the detector with the dark matter galactic wind changes; so that phonons with different momenta—i.e.aligned to the dark matter wind orientation—are preferentially generated during different times of the day. Moreover, as the DM interaction is bound to the crystalline symmetry through the collective excitations in the solid (phonons), two nonequivalent crystalline directions will interact differently with the DM. Consequently, for anisotropic crystals, a daily modulation of the signal is expected, due to the different interaction strengths for the wind-oriented initial states. This particular feature has been proposed as a way to discern a DM generated signal from other sources and/or noise; nevertheless the effect and possible masking by downconversion processes of such an interesting feature for detection has never been discussed.

In Supplementary Figures 2-3 we present the signal for initial states with momentum parallel and perpendicular to the \(\langle0001\rangle\) direction, which correspond respectively to the ordinary and extraordinary DM interaction axis of Al\(_2\)O\(_3\) [46]. It is clear from those, that phonon downconversion together with interface scattering yield similar fluxes for both—that is perpendicular and parallel—initial states. Therefore, any initial imbalance in the phonon generation rate would directly translate to the collected flux. In other words, the daily modulation in the Al₂O₃-\(\langle111\rangle\) Al scheme is not masked either by phonon dynamics within the target or by the target-collector interface scattering.

Feasibility of phonon-based detection of sub-MeV DM. Up to this point we have studied the influence of various parameters –the phonon energy, the phonon momentum, the presence of isotopes, where the athermal phonon is generated within the target– on the signal that reaches the collector and we have obtained a set of design rules for an efficient detector. Now we assess the overall feasibility of the detection scheme considered. In other words, can a carefully optimized detector provide a sufficiently intense signal that a TES can detect?

To this end, we have extended one of our Monte Carlo simulations up to 1 \(\mu\)s, focusing on the Al₃O₂\(\langle0001\rangle\)/Al\(\langle111\rangle\) case with \(L=\SI{10}{\nano\meter}\), \(E=\SI{0.06}{\electronvolt}\), \(\lVert\mathbf{q}\rVert=\SI{1}{eV}\), and \(d=\SI{2.5}{\nano\meter}\). Then we plotted the power as a function of time and compared to the sensitivity of state-of-the-art TESs (we have assumed \(\mathrm{NEP} = 1.5\times 10^{-18}\,\mathrm{W}/\sqrt{\mathrm{Hz}}\) and \(\mathrm{bandwidth} = 2.6 \,\mathrm{kHz}\)) [32], [71]. The results are shown in Figure 11. As mentioned above, TES have rather slow response times, \(t_r \gtrsim \mu\mathrm{s}\), and this plot highlights how critical this parameter is. The signal would be in the limit of detectability for state-of-the-art single TES devices  [32], [72], indicating the need for further improvement of both time response and NEP of TES. The detection scheme, though, seems feasible either with the progress of present sensors or with new designs [73].

A related question is whether sufficient target volume can be achieved to be sensitive to well-motivated models of Dark Matter. Many studies have provided a simple estimate of the sensitivity by calculating the exposure which is required in order to produce \(\mathcal{O}(1)\) DM-induced phonons [8], [15], [33], [39], [47]–[49], [74], [75]. Calculating the DM-phonon scattering rate using  [46], we find that currently unconstrained parameter space can be explore with an exposure \(\sim 10^{-3}\,\text{g-days}\), while models in which DM is produced with the correct relic abundance via freeze-in would require an exposure of \(\sim 1\,\text{g-days}\) [15]. Considering a few days of exposure and a surface area of \(\sim 100\,\mathrm{cm}^2\) for the target crystal, this would require target thicknesses of \(\mathcal{O}(10\,\mathrm{nm} - 10\,\mu\mathrm{m})\). Of course, a detailed calculation is required, including the phonon energy losses described in this work, but this estimate demonstrates that wafers of realistic thicknesses can be a feasible target for probing unexplored DM parameter space.

Finally, we find it important to highlight that our results should be understood as an upper limit, as a real device will contain further sources of signal diminution. These include sample defects or the partial coverage of the target by the collector that has not been considered in our simulations, not to forget about losses related to other processes within the collector that are not related with phonon dynamics. Also, our simulation scheme contains several approximations to reduce the computational burden (e.g. use of bulk properties, elastic scattering at the interface, etc.) that can further affect/diminish the signal (see Appendix 4.3 for a comprehensive list and discussion of the approximations and limitations of the applied methodology).

4 Conclusions↩︎

We have presented a detailed scrutiny of the phonon dynamics resulting from an athermal optical phonon excitation, followed by downconversion and propagation within a semiconductor target, until reaching the interface of a phonon collector. We consider these processes in the context of the detection of light DM, though they are general to a wider class of instruments based on cryogenic thermal detectors. Our formalism is entirely based on parameter-free, ab initio electronic structure calculations. We focused on an Al\(_2\)O\(_3\)/Al target/collector, two material systems that offer several advantages for the detection of light DM, such as a phonon spectrum of the target in the right energy range, an anisotropic response suited to reveal daily modulation and a high-quality interface between the semiconductor and the superconductor.

A preliminary methodological survey revealed that any RTA-based approach is bound to yield predictions of limited accuracy, largely underestimating the heat flux that reaches the collector and then artificially dictating experimental conditions that are more difficult to meet. Conversely, linearization of the PBTE works well. Based on these conclusions, we developed a 3D beyond-RTA phonon Monte Carlo that allowed us to introduce the spatial dimension of the device and address questions about its size or the effect of where in the target the DM scattering takes place (i.e., the distance from the collector). We showed that the wavevector of the optical phonon initially excited by the DM particle plays a negligible role, while a higher energy phonon results in a (slightly) larger heat flux reaching the collector. Isotopes can, perhaps counterintuitively, favor the detection and result in a larger heat flux by providing transport channels of higher velocities. Analyzing these effects in materials with a larger dispersion in the isotope population (e.g., SiC) deserves a separate study to see if a more conventional behavior is recovered (i.e., isotopes decrease the heat flux) or if the behavior reported here is amplified. The mismatch in the phonon spectra of Al\(_2\)O\(_3\) and Al is the source of a considerable thermal boundary resistance, which nonetheless does not limit significantly the efficiency of the detector: phonons undergo multiple reflections between the target/collector interface and the target backsurface, until they fit the narrower energy window of the collector. For this reason, too thick targets should be avoided (the multiple reflections would imply too large distances travelled, with the consequent losses), while the initial position of the phonon excitation does not seem to play a major role.

Altogether, our results suggest that, though challenging—particularly in relation to the sensitivity of the TES with respect to the collected power—the direct detection of light DM via athermal phonon generation using a TES-based detector appears to be feasible, or at least, the phonon downconversion followed by quasi-ballistic phonon propagation does not appear to be a major bottleneck in terms of reducing the signal. We also stress that most of these conclusions are general and can be extended to any phonon-based cryogenic detection scheme, regardless of the source (DM in our case) of the initial phonon excitation. On the downside, many of the results presented here are expected to be material-dependent and thus the quantitative importance of the different factors should be carefully evaluated case by case.

ACKNOWLEDGMENTS↩︎

We acknowledge financial support from AEI under grant PID2020-119777GB-I00, and the Severo Ochoa Centres of Excellence Program under grant CEX2019-000917-S, by the Generalitat de Catalunya under grant no. and 2017 SGR 1506, and by the Spanish NRRP through the CSIC PTI QTEP+. BJK acknowledges funding from the Ramón y Cajal Grant RYC2021-034757-I, financed by MCIN/AEI/10.13039/501100011033 and by the European Union “NextGenerationEU"/PRTR. Calculations were performed at the Centro de Supercomputación de Galicia (CESGA) within action FI-2022-3-0031,”Phonon-based detection of ultralight dark matter”, of the Red Española de Supercomputación (RES). The authors also thank the ‘Dark Collaboration at IFCA’ working group for useful discussions.

APPENDIX: METHODS EXTENDED↩︎

4.1 Solving the Swartz’s model↩︎

The key ingredients of Swartz’s model (Eq. 5 ) are the time evolution of the phonon deviational population [\(n^d(t)\)], and the transmission probability (\(\alpha\)) The former, under the assumption of homogeneity—as it is the case of Swartz’s model—, can be obtained by solely time integrating the collision operator:

and/or its variations, i.e. the linearized scattering operator—that is obtained by neglecting the order-2 deviational terms—and the pseudo-RTA (see Eq. 4 ). We note that Eq. 8 only includes isotopic and three-phonon scattering process.

Generally, this equation contains two or more time scales of a different order, i.e. some phonon modes are evolving much faster than others, a characteristic typical of a stiff system. So that, the optimal methodology to integrate Eq. 8 would be to make use of an implicit method (e.g. Rosenbrock methods or implicit Runge-Kutta methods); nonetheless, such methods require the computation of the Jacobian. However, the computation of such a matrix is computationally prohibitive in our case; for reference, in cubic GaAs, taking a \(24\times24\times24\) \(\mathbf{q}\)-mesh would mean a Jacobian of \(\sim\)​55 GB. Consequently, one is constricted to make use of explicit methods with a small enough timestep (\(\Delta t\))—either arbitrarily selected by the user or chosen based on the RTA lifetimes, namely \(\Delta t = \min\{0.25\mathrm{E}(\tau_i|\tau_i \neq 0),1\}\) ps—to integrate \(n_i(t)\).

Amid all the available explicit methods, the predictor-corrector methods are especially suited for problems in which the computation of the derivative is computationally expensive, as it is the case of the full collision operator. Therefore, we have selected a two-step predictor-corrector method with a constant timestep as implemented in boost::odeint [76] library, to integrate Eq. 8 and its variations, using an explicit \(4^{\mathrm{th}}\) order Runge-Kutta to obtain the initial states required for the selected scheme. We note, that for the pseudo-RTA case it is much more efficient to directly solve the Eq. 3 ; so that the deviational phonon population evolves accordingly to \[n^d(t+\Delta t) = \exp(A \Delta t) n^d(t) \approx (I + A \Delta t) n^d(t),\] where we have expanded the exponential matrix [\(\exp(A \Delta t)\)] using a Taylor series up to the 1^st order.

In the building of this collision operator it is essential for it to respect crystal symmetries and microscopic reversibility [56]; consequently we have used the symmetric adaptive smearing scheme of Ref. [77] to compute the three-phonon processes.

Besides symmetry, it must also be noted that the linearized and the full scattering operator must in principle conserve the energy (note again here that RTA, and pseudo-RTA are energy destructive by construction); however, this condition is violated due to the broadening scheme, in which delta conservation energies are approximated by Gaussians. Therefore, we add a correction (\(\beta_i\)) extracted from a Lagrange-multiplier approach, similarly to the one developed in Ref. [56] to our derivative computation each step to enforce energy conservation. So that, given a violation of energy, \(\Delta E = \sum_i \hbar \omega_i\frac{\partial n^d_i}{\partial t}\), imposes a constraint for the applied correction \(\left( \sum_i \beta_i = -\Delta E \right)\). Moreover, in order to reduce possible side effects of such a correction, \(\beta_i\) must remain as small as possible. This can be achieved through a Lagrange multiplier’s through the optimization of the following constrained objective function: \[F = \sum_i \beta_i^2 + \lambda\left(\sum_i \beta_i + \Delta E\right),\] and thus, \[\begin{gather} \frac{\partial F}{\partial \beta_j} = 2\beta_j + \lambda , \end{gather}\] which is minimized by \(\beta_j = -\frac{\lambda}{2}\). Consequently, combining this last with the constraint, one arrives to: \[\beta_i = -\frac{\sum_j \hbar \omega_j\frac{\partial n^d_j}{\partial t}}{N_{\mathrm{modes}}}\] where \({N_{\mathrm{modes}}}\) is the number of modes.

Finally, as for the latter ingredient, the transmission coefficient, we compute it using as previously mentioned the DMM; nevertheless we also allow for a uniform setting of an equal transmission coefficient for all the modes. This last is a very simplistic model that works for pairs with similar maximum frequencies like GaAs-Si, but yields unreasonable results for most cases, like the present case, as it would allow transmission from high-energy modes of Al\(_2\)O\(_3\) that do not exist in Al.

4.2 Practical notes on MC implementation: the propagator↩︎

The basics of the MC solver can be found on Ref. [77]—BTE-Barna, which this code is derived from—, and more theoreticall background can be found in Refs. [56] and [78]. Here we only highlight the main difference with respect to BTE-Barna, that is the computation of the propagator, namely \(\exp(B\delta t)\), required for the scattering algorithm [78]. In contrast with previous works [56], [77], [78], we are dealing with 3D materials instead of 2D ones, rocketing the memory and computational resources required for the computation; therefore, it becomes desirable to keep up not only the sparsity of the transition matrix but in the propagator itself.

Consequently, we introduce the energy conservation correction only on non-zero entries, so that the transition matrix sparsity is conserved, but we also make use of the first-order Taylor expansion to the exponential—i.e. \(\exp(B\delta t) \approx I + B\delta t\) to translate this sparsity to the propagator itself. Nonetheless, this last approximation can seem crude, it was proven by Landon that for small time steps—namely like those used to solve the Swartz’s model—is virtually equivalent to more refined methods [78].

4.3 Practical notes on MC implementation: approximations and limitations↩︎

Although our approach is in principle the most exact one among those used up to now, our simulations have still several limits and approximations that should be kept on mind, most of them coming from a design point of view. Here we highlight some of these limitations.

• Boundary scattering: The implementation is like the one of the RTA MC solvers included in BTE-Barna, i.e. perfectly elastic and diffusive.

• Spatial resolution on in-plane direction: The in-plane spatial discretitation is disregarded, and the system is assumed to be infinite in those directions. This last can be only considered if the transport direction (out-of-plane) is much smaller than the others.

• Phonons in the collector: we are assuming that phonons do not reenter from the collector to the target. This last is partially justified by them being absorbed in breaking the Cooper pairs in the collector.

• Phonon properties in the target/collector: we are assuming bulk phonon properties, both harmonic and anharmonic in both target and collector.

References↩︎