Microwave Photon Number Resolving Detector Using the Topological Surface State of Superconducting Cadmium Arsenide

Photon number resolving detectors play a central role in quantum optics. A key challenge in resolving the number of absorbed photons in the microwave frequency range is finding a suitable material that provides not only an appropriate band structure for absorbing low-energy photons but also a means of detecting a discrete photoelectron excitation. To this end, we propose to measure the temperature gain after absorbing a photon using superconducting cadmium arsenide (Cd3As2) with a topological semimetallic surface state as the detector. The surface electrons absorb the incoming photons and then transfer the excess energy via heat to the superconducting bulk's phonon modes. The temperature gain can be determined by measuring the change in the zero-bias bulk resistivity, which does not significantly affect the lattice dynamics. Moreover, the obtained temperature gain scales discretely with the number of absorbed photons, enabling a photon-number resolving function. Here, we will calculate the temperature increase as a function of the number and frequency of photons absorbed. We will also derive the timescale for the heat transfer process from the surface electrons to the bulk phonons. We will specifically show that the transfer processes are fast enough to ignore heat dissipation loss.

Photon number resolving detectors play a central role in quantum optics. A key challenge in resolving the number of absorbed photons in the microwave frequency range is finding a suitable material that provides not only an appropriate band structure for absorbing low-energy photons but also a means of detecting a discrete photoelectron excitation. To this end, we propose to measure the temperature gain after absorbing a photon using superconducting cadmium arsenide (Cd3As2) with a topological semimetallic surface state as the detector. The surface electrons absorb the incoming photons and then transfer the excess energy via heat to the superconducting bulk's phonon modes. The temperature gain can be determined by measuring the change in the zero-bias bulk resistivity, which does not significantly affect the lattice dynamics. Moreover, the obtained temperature gain scales discretely with the number of absorbed photons, enabling a photon-number resolving function. Here, we will calculate the temperature increase as a function of the number and frequency of photons absorbed. We will also derive the timescale for the heat transfer process from the surface electrons to the bulk phonons. We will specifically show that the transfer processes are fast enough to ignore heat dissipation loss.

I. INTRODUCTION
Photon number resolving detectors have been explored significantly over the past decades [1][2][3] due to the dire need for resolving the number of photons in applications such as the security of quantum communications [4,5] and the sensitivity of quantum sensing [6]. As photon-based quantum computing advances, precise resolution of photon number detection is increasingly important. Microwave photons are the backbone of prolific transmon quantum computation, and therefore, detection of the microwave photons is tremendously important in the current quantum computing paradigm [7]. Several non-number-resolving techniques to detect microwave photons have been developed, including the circuit QED technique [8], dressed-state superconducting quantum circuit [9], current-biased Josephson junction [10], and the dark-state detector [11]. It is well-known that building a parallel detection system, first splitting the light path using beam splitters and then using nonnumber-resolving detectors in each parallel path, may provide a probabilistic photon number resolving detection, which is further limited due to the loss associated with parallelization. In contrast, a single photon-number resolving detector with a deterministic photon number resolution would provide a immense advantage particularly in photonic quantum computers by reducing the error-correcting overhead. To the best of our knowledge, a single-device photon-number resolving detector that can simultaneously detect multiple incoming photons at microwave frequency has never been reported so far.
Here, we propose a photon-number resolving detector operating at microwave frequencies, based on the topological surface states of cadmium arsenide (Cd 3 As 2 ). Semimetals such as graphene provide an ideal detecting material for microwave photons due to their zero band gap. Recently, Dirac and Weyl semimetals with Dirac cone dispersion have gained prominence due to high mobility [12], along with the fact that they can be synthesized through conventional techniques [13][14][15]. Particularly, Cd 3 As 2 displays proximity-induced bulk superconductivity at low temperatures, and the electronic structures of the bulk and the topological surface states are decoupled. Maintaining the Cd 3 As 2 semimetal material at a very low temperature is necessary for an efficient photon-induced electron excitation to a conduction band just above the Fermi level due to the low photon energy. The bulk state enters a superconducting state at a sufficiently low temperature, opening a band gap beyond the microwave photon energy. Fortunately, the topological surface state of Cd 3 As 2 is not affected by the temperature, continuining to provide a gapless Dirac cone. We use this topological surface state as a photon absorber. Once the photon is absorbed, a rapid rethermalization in band population occurs with a new elevated temperature corresponding to the absorbed photon energy. We then utilize the fact that the redistributed electron population transfers its energy to the bulk's phonon modes via a surface electron-bulk phonon coupling, thus increasing the bulk's temperature. The elevated bulk temperature then reduces the conductance of the superconducting bulk electron state, which is measured and used to eventually indicate the number of photons absorbed.
The paper is organized as follows. In Sec. II, we briefly review the photon absorption in the topological surface state of Cd 3 As 2 . In Sec. III, we connect the event of photon absorption in the topological surface state electrons to the bulk temperature increase via a two-step process, namely, the energy gain for surface electron modes and, then, the transfer of electron energy to the bulk phonon modes. Sec. IV resolves the important issue of the time scale of the temperature increase and shows that the ab-FIG. 1: Basic layout for the Cd 3 As 2 photon number resolving device. At low temperatures, the bulk becomes superconducting, while the (112) surface retains a graphene-like dispersion. A photon (depicted by the green arrow) is absorbed by the surface electrons. The change in bulk resistivity (measured by the electrodes at zero bias) is used to determine the temperature increase. sorbed photon energy indeed is most likely transferred to the increase of the bulk temperature, rather than being lost to radiative decay of the excited electrons. Sec. V presents the construction of a photon-number resolving detector based on the results of the previous sections. Numerical examples with realistic parameters build credible real use cases of the proposed scheme. In the final section, we summarize our results and suggest the path towards building a photon-number resolving detector with near-unity efficiency on a chip.

II. PHOTON ABSORPTION IN TOPOLOGICAL SURFACE STATE
The setup for the device is depicted in Figure 1. The system is based on a Cd 3 As 2 crystal inside a low temperature refrigerator with a baseline temperature below the bulk superconducting critical temperature. Recent experimental findings have demonstrated that Cd 3 As 2 features a topological surface state on the (112) surface with a linear band crossing around the Dirac points [16,17]. This graphene-like surface state band structure at low energy can be attributed to the fact that sites consisting of stacked As and Cd atoms approximately form a honeycomb superlattice on the (112) surface [13]. As with graphene, the dispersion relationship can be expressed as the following linear function of the wavevector k when the Fermi level is at the Dirac point: where is the Planck's constant, c and v represent the conduction and valence bands, respectively, and v F denotes the Fermi velocity, which is approximately 10 6 m/s [16,18] its semimetallic property [19]. The superconducting gap frequency f ∆ corresponding to this critical temperature T c can be determined using the standard BCS equation [20]: where k B is the Boltzmann constant. The resulting band structures for the superconducting bulk and the topological surface state are compared in Figure 2. Photons of microwave frequency below 50 GHz will be absorbed solely by the surface state, as desired.
In order to derive the absorption probability, we consider the physical manifestation of the Dirac cone on the nature of the Bloch states. As in graphene, each electronic state in the vicinity of a Dirac point can be conceptualized as a massless Dirac fermion with a well-defined momentum k, where k represents the wavevector of the electronic state in the reciprocal space for which the Dirac point is the origin. Therefore, the absorption coefficient for Cd 3 As 2 will equal the well-known corresponding value for graphene [21]: where e is the charge of an electron, 0 is the vacuum permittivity, c the speed of light, and f (E) denotes the Fermi-Dirac electron occupation probability at energy E. Note that the absorption coefficient (as a function of photon frequency) is invariant with respect to the Fermi velocity. This is because the interband dipole matrix element (corresponding to the absorption probability for a single electronic state) increases with the Fermi velocity, while the density of states decreases with it, eventually canceling each other's effect. At high frequencies, the absorption coefficient is approximately invariant with frequency, equaling a constant value of e 2 /(4 0 c) = 2.3%. On the other hand, at low frequencies on the scale ω k B T , the coefficient becomes attenuated, reaching a minimum value of one-fourth the high-frequency absorption as the Dirac cone is approached.

III. TEMPERATURE INCREASE VS. ABSORBED PHOTON NUMBER
Having determined the probability that the topological surface state absorbs a photon from an incoming field, our next step is to determine how the absorption of a single photon increases the temperature of the 3D bulk sample. When a photon excites an electron to the conduction band, a rapid rethermalization of the Fermi-Dirac distribution through electron-electron interaction ensues, leading to an electron temperature above the lattice temperature. For undoped monolayer graphene, which features a band structure approximately identical to that for the Cd 3 As 2 surface state, this process occurs in tens of picoseconds for cryogenic baseline temperatures [22]. The carrier rethermalization is followed by heat transfer from the collection of electrons to the lattice via electron-acoustic phonon interaction, until a thermal equilibrium is reached between the electron temperature and the lattice temperature. Generally, the interaction between electrons and acoustic phonons is much slower than the electron-electron interaction [23][24][25]. Afterward, the Fermi-Dirac distribution for the electron bands and the Bose-Einstein distribution for the phonon branches will comply with the same temperature.
We now determine the temperature increase due to the absorption of a photon of frequency ω by calculating the portion of the imparted energy that is eventually converted to bulk lattice vibrations (i.e. the phonon modes) and to the surface electron modes, and by deriving the heat capacity of the two systems. We will make two important assumptions here: firstly, that the electronelectron interaction rate dominates over the radiative loss rate of the electrons (which we will demonstrate in a later section), and secondly, that the very low values for the initial and final temperatures ensure that virtually all of the phonons are located in the low-energy linear parts of the acoustic branches, thus allowing for use of the Debye approximation [26] in determining the heat capacity.

A. Energy Gain for Surface Electron Modes
We start by writing out an expression corresponding to energy conservation in the system given the absorption of a photon of frequency ω: where ∆U el and ∆U ph represent the total energy gained by the topological surface electron modes and the bulk phonon modes, respectively, at equilibrium. We focus first on the energy gain for the electron modes as a function of electron temperature, as this will be necessary for calculating the initial electron temperature gain after photon absorption but prior to heat transfer to the lattice. The total electron energy with respect to the Fermi sea is calculated by taking a sum of the conduction and valence band energies weighted by the Fermi-Dirac occupation probabilities (multiplied by 2 to account for the fact that each spatial state contains 2 spin states): (5) The first term corresponds to the energy gained in creating a conduction band electron, while the second term corresponds to the energy gained in creating a valence band hole. Since the deviation from the Fermi sea at low temperatures is concentrated in the vicinity of the Dirac cone, we can assume that the linear isotropic dispersion relationship in Eq. (1) holds for the entire relevant wavevector range. Therefore, the summation over wavevectors can be replaced by an integral over the density of spatial states (ρ c and ρ v for the conduction and valence bands, respectively): Since the conduction and valence band energies are opposite at each wavevector, we can express U el as a single integral over the energy absolute value E, where E c = E and E v = −E. Due to the equivalent magnitudes of the dispersion slope for the conduction and valence bands at each wavevector, we can further define a general density of states The density of states at band energy E can be solved by applying the dispersion relationship as follows: where A is the surface state area, A k is the area in the reciprocal space associated with the value k = |k|, and N is the number of states. As expected for the graphenelike band structure, the density of states is linear in the energy. We are thus in a position to solve for the electron energy U el (T ) as a function of temperature T : where ζ represents the Riemann zeta function, with ζ(3) ≈ 1.2. This implies that the electron temperature varies with the total electron energy as U 1/3 el , and the relationship between the gain in electron energy and the temperature change from T i to T f takes the following form: For detecting a moderate number of microwave photons, we are primarily interested in the limit ∆T In that limit, the cooling power Q for the electron modes is related to the rate of change of the temperature by taking the derivative of ∆U el with respect to time, yielding the following function of the temperature T ≈ T i :

B. Energy Gain for Bulk Phonon Modes
Next, we look to determine how the energy gained by the lattice vibrations relates to the lattice temperature. In general, the total phonon energy is determined as a function of temperature by summing over the modes corresponding to various phonon branches µ and phonon wavevectors q, weighted by the occupation number n µ,q for each phonon mode: where ω µ,q is the frequency of the phonon mode of branch µ and wavevector q, and the occupation number at a given temperature T is calculated from the Bose-Einstein distribution: As previously mentioned, the fact that the sample is in the low temperature regime (below 0.5 K) indicates that the Debye model, with its assumption of a linear phonon dispersion, is approximately valid for the phonon modes with non-negligible occupation numbers. Therefore, we can restrict the summation over branches to just the 3 acoustic branches (corresponding to the 3 polarizations).
In general, the slope of each of these branches is slightly anisotropic in reciprocal space due to the varying angle between the polarization and propagation directions. However, per the treatment in Kittel [27], we can approximate the composite effect of the 3 branches on the summation in Eq. (12) as equivalent to a summation over 3 isotropic branches, each featuring a slope of v s (i.e. the speed of sound in the material), such that: where ω q = v s q, and the speed of sound v s for Cd 3 As 2 is estimated as 2.3 × 10 3 m/s [25]. We now replace the summation over wavevectors with an integral over the density of states for each branch in terms of frequency, D(ω). For a 3-dimensional lattice with a speed of sound v s , this density is determined as follows: where V is the bulk volume of the lattice and V q is the volume in the reciprocal space associated with q = |q|. We therefore solve for the total phonon energy as a function of temperature through the following integral: Note that the expression differs from that in Kittel in that we use ∞ instead of a specific Debye cutoff for the upper bound of the frequency range. As in the case of the integral over energy for the electronic modes in the previous section, this is justified by the rapid convergence of the integrand to 0 at very low temperatures [28], corresponding to the fact that only the linear regime of the acoustic branches are non-negligibly occupied. Then, we find the following result: Unlike the total surface electron energy, which scales as T 3 , the total bulk phonon energy scales as T 4 . This difference can be attributed to the fact that the surface is 2D whereas the bulk is 3D, having more degrees of freedom, thus implying that all else being equal, a given change in energy would have a weaker effect on bulk temperature than on surface temperature. The change in the total phonon energy can therefore be related to the initial and final temperatures T i and T f as follows: Note that, for a prism-shaped sample (such as a thin film) for which the surface state forms one of the two bases for the prism, the bulk volume is proportional to the surface area as V = Ad, where d represents the sample thickness.

C. Comparison of Specific Heat
Having derived the energy gain for the surface electron and the bulk phonon modes for given initial and final temperatures, we now seek to compare the specific heat values for the two mode types in order to glean an understanding of how the excess thermal energy is distributed between the modes. From the result in Eq. (10), the specific heat for the collection of surface electrons is determined as follows: The specific heat for the bulk lattice is calculated in an analogous manner from the result in Eq. (18): Using the relationship V = Ad, where d denotes the bulk depth of the lattice, we divide Eq. (20) by (19) in order to determine the ratio of the specific heat values for a given temperature: Since the baseline refrigerator temperature is at least 0.25 K, this sets the minimum value for T . The lattice depth d must be at least multiple times longer than the lattice constant, which is 3 − 5Å for Cd 3 As 2 [29]. Therefore, the phonon specific heat is far higher than the electron specific heat (by at least 3 orders of magnitude), indicating that nearly all of the thermal energy gained from the photon absorption is eventually stored in the lattice vibrational modes. As such, if N photons of frequency ω are absorbed, then the relationship between the final equilibrium temperature T f and the initial temperature T i can be determined from Eq. (18): As this expression shows, the determinitive factor in the temperature increase is the total photon energy absorbed per unit volume of the lattice, which is proportional to N ω/V . We plot the temperature gain as a function of photoelectron density N/V for 3 different baseline temperatures 0.35 K, 0.40 K, and 0.45 K, for an input photon frequency of 5 GHz (corresponding to ω = π × 10 10 s −1 ) in Figure 3. Note that for low photon densities such that the temperature gain is small compared to the baseline temperature, the relationship between temperature gain and photon density is approximately linear, as expected.
The imbalance between the electron and lattice specific heat values also has significant implications for the heat transfer between the electron and phonon distributions that ultimately yields the equilibrium state. In particular, when comparing the quasiequilibrium electron and lattice temperatures to the final equilibrium temperature for the whole system, the equilibrium temperature will be much closer to the quasiequilbrium lattice temperature than to the electron temperature. We consider the timescale for the electron-lattice heat transfer as a func-tion of temperature in the next section.

IV. ELECTRON-PHONON INTERACTION TIMESCALE
Having derived the equilibrium temperature for the bulk lattice upon photon absorption, we now aim to estimate the timescale over which that equilibrium is reached. As previously mentioned, this energy transfer takes place in two steps: a rapid electron-electron rethermalization, followed by heat transfer from the electrons to the acoustic vibrations of the lattice (which is much slower than the electron-electron interaction [23][24][25]). Here, we will focus on the latter process, since it serves as the limiting factor in setting the minimum timescale for reaching equilibrium. We will characterize the available phase space area for bulk phonon emission by the surface electrons, derive the matrix element for the electron-phonon interaction, and finally calculate the rate for the electron-phonon heat transfer.

A. Phase Space
We start by examining the available phase space for the interaction between 2D surface state electrons and the bulk phonon modes. Unlike the spherical equal-energy manifolds for the 3D Dirac cone carrier modes, the 2D Dirac cone electron modes take a cylindrical equal-energy manifolds, with a degree of freedom in the k z -direction (corresponding to the axis perpendicular to the surface). We therefore use cylindrical coordinates, expanding the final electron wavevector p as (p cos θ p , p sin θ p , p z ) and the initial wavevector k as (k, 0, 0). In this coordinate system the emitted phonon wavevector q can be expanded as follows: Note that the bulk phonon's equal-energy manifolds retain a 3D spherical shape defined by ω = v s q. We map this onto the electron's manifolds using energy conservation: Squaring both sides and solving for p, we find that p varies with both θ p and p z : where ∆(θ p ) is defined as follows: For our material, v s v F , thus yielding ∆(θ p ) 1 for all θ p . Therefore, p can be approximately expressed as solely a linear function of p z : Geometrically, the available phase space can be envisioned as pair of cones aligned along the p z -axis, with the bases overlapping at p z = 0, as depicted in Figure 4. The radius attains is maximal value of k at p z = 0 and tapers off as the magnitude of p z increases. The vertices are reached at the following values of p z : Finally, we determine the phonon wavevector q from the calculated value of p as a function of k. As shown in Eq. (23), the amplitudes of p z and q z will equal each other, since the electron and phonon dispersion centers lie on the same xy-plane. We label q xy as the component of q perpendicular to the q z -axis. For a given p z and θ p pair, the amplitudes q xy , k, and p form the 3 legs of a triangle for which θ p represents the angle between the sides of lengths k and p. Therefore, q xy can be calculated as follows: and since q z = −p z , the frequency ω q of the emitted phonon is straightforwardly calculated from the speed of sound: Since v F v s (by a factor of 400), the length of each cone is far longer than the diameter, implying that the approximation q ≈ |q z | = |p z | ≈ p will be valid for nearly all of the available phase space.

B. Heat Transfer Rate
Having determined the phase space for the electronphonon interaction, we are now in a position to calculate the heat transfer rate between the two modes from the composite interaction. Labeling the electron energy for a generic wavevector k as E k , the matrix element corresponding to the electronic transition from k to p through the emission of a phonon in branch µ and wavevector q as M µ,q k,p , the Bose-Einstein phonon occupation number for the mode frequency ω at temperature T as n T (ω), and the Fermi-Dirac distribution value at T as f (T ), the rate Q is determined through the following summation over initial carrier wavevectors k, final carrier wavevectors p,  = (k, 0, 0). The phase space forms a double cone, with a base of radius k on the p x p y -plane, and tapering off in the +p z and −p z directions with a length of v F k/v s along each. Note that the base vs. height ratio for the double-cone in (b) is not to scale. and phonon branches and wavevectors (µ, q) [30,31]: As previously discussed, the low temperature restricts the occupied phonon modes to the long-wavelength acoustic regime. The interaction between electrons and long-wavelength acoustic phonons is dominated by the deformation potential [32,33], as recently applied to the interaction between bulk electrons and phonons in Cd 3 As 2 [25,34]. As discussed in Appendix A, the equivalent matrix elements apply for the interaction between surface electrons and bulk phonons. Therefore, the matrix element amplitude-squared reduces to a function varying solely with and linear in the phonon amplitude q: where V represents the lattice volume and C is a constant that varies with the square of the deformation potential. Substituting this, along with the electron and phonon dispersion relationships into the expression for Q, we find that it takes the following form: The summation is simplified in Appendix B in the limit ∆T T , where T ≈ T e ≈ T L and ∆T is defined as T e − T L . We find that the Dirac and Kronecker delta functions combine to reduce the integral over the phase space volume to the double-cone phase space area derived previously, as expected. This yields the following expression for the carrier-phonon heat transfer rate due to intraband (valence-valence or conduction-conduction) transitions: Solving the integral numerically, we obtain a value of −32. Therefore, Q is further reduced to the following: Next, we solve for the heat transfer rate due to interband transitions. Based again on Appendix B, we use the following expression: Note that the constants in front of the integral are identical to that for the intraband case. Solving this integral numerically yields a value of −100. The total heat transfer rate Q total from the surface carriers to the lattice vibrations is determined by multipyling the intraband rate Q by 2 (to account for both bands) and then summing with the interband rate Q inter : In order to determine the heat transfer timescale, we substitute the previously derived relationship between the electron cooling rate and the rate of change of electron temperature from Eq. (11) into the left-hand-side of the above expression: As the result shows, the electron temperature decays exponentially toward the lattice temperature, with the rate varying as T 3 . The remaining task is to determine the value of the constant C, which derives from the electron-phonon matrix element. One method for doing so is by using the deformation potential of 20 eV measured by Jay-Gerin et al. [35]. This yields the following value for C, using a material density ρ = 7 × 10 3 kg/m 3 [25]: This leads to the following heat transfer time constant γ: An alternative method for finding the deformation potential is by merging the experimental results from Weber et al. [36] with the theory provided by Lundgren and Fiete [25]. Specifically, Weber et al. used a bulk Cd 3 As 2 sample intrinsically doped to a baseline electron density of 6 × 10 23 m −3 , which corresponds to a Fermi energy of 170 meV and a Fermi temperature of 1130 K. Under these conditions, they observed a timescale of 3.1 ps for electron cooling by low-energy acoustic phonon emission at lattice temperatures of 80 K and 300 K. This scenario is addressed by Lundgren and Fiete's Equation (8), which models the heat transfer rate for k B T E f (where E f is the Fermi energy): We now substitute a temperature and rate data point from Weber et al. into this expression to calculate the deformation potential D. Since the limit k B T E f is much more valid for T = 80 K than for 300 K, we use the former as the temperature corresponding to the rate γ = (3.1 ps) −1 for the purposes of application to Eq. (41). This yields the following value for D: This leads to the following value for the coefficient C: which yields the following heat transfer time for our model: As will be discussed in the next section, the lower bound for the baseline temperature T (which will also set the minimum value for the heat transfer rate) will be about 0.35 K. For this temperature, the above two methods yield a lattice heating timescale approximately ranging from 93 ns to 15 µs.

V. PHOTON-NUMBER RESOLVING DETECTION
We now describe the photon-number resolving detector scheme based on our theoretical findings. First, we address the question whether the timescale for lattice temperature equilibration is much faster than the dissipation time through thermal conduction or radiative decay. Regarding the thermal conduction heat loss, we note that the contacts used for cooling the sample can be removed after the material reaches the refrigerator temperature. As a result, the heat dissipation time through thermal conduction will range on the order of several hours and can thus be ignored. Instead, we focus on the radiative loss. Based on the results calculated for graphene, the electron-hole interband dipole moment for a 2D Dirac cone band structure is given as a function of photon radial frequency ω as follows [37]: Substituting this into the well-known radiative decay rate expression based on the Einstein coefficients [38], we find that the radiative rate varies linearly with ω: For frequencies up to 10 GHz, the radiative decay time is therefore 150 µs or greater. This is significantly longer than the electron-phonon heat transfer time calculated above, which is 15 µs or less, which in turn is much longer than the previously discussed electron-electron rethermalization time of tens of picoseconds [22]. Therefore, a rapid rethermalization of the electron population in the bands occurs before any radiative loss of the photoelectrons occurs. We thus conclude that the energy of the absorbed photons is safely transferred to, firstly, the rethermalization of the carrier band populations, and then, to the bulk phonon modes to elevate the bulk temperature. Next, we discuss how the bulk temperature is measured. Since the elevated bulk temperature will increase the bulk resistance of the superconducting bulk states as shown in Figure 5, we measure the zero-bias resistivity across the bulk (using a lock-in amplifier) as a proxy for the temperature. This is advantageous relative to infrared-based bolometry since it does not perturb the electronic structure of the bulk, as well as due to the fact that electrical signals can be measured in ultrafast picosecond-range intervals [39]. We manufactured a Cd 3 As 2 device to measure the superconducting bulk resistivity as a function of sample temperature. To this end, it is important to note the lower bounds for the dimensions of each Cd 3 As 2 crystal. The goal of the device is to measure photons in the transmon frequency range, i.e. 5 − 7 GHz [40,41]. For a Dirac cone dispersion with Fermi velocity v F , a photon of frequency f is resonant with the band gap at the following band wavevector: Therefore, in order for resonance to exist at photon frequencies as low as 5 GHz, the maximum length of each Bloch state in reciprocal space must be ∆k ≈ 1.6 × 10 4 m −1 , thus implying that the minimum length of the Cd 3 As 2 surface along each dimension is 2π/∆k = 0.4 mm. We also assume that the depth of the lattice is limited by design constraints to a minimum value of 20 nm, since this is the minimum thickness that has been achieved with an MBE technique [14]. For a photon frequency of 5 GHz and crystal dimensions of 0.4 mm by 0.4 mm by 20 nm, the single-photon temperature gain is calculated by substituting the values N = 1, ω = π × 10 10 s −1 , and V = 3.2 × 10 −15 m 3 into Eq. (22) and linearizing: For temperatures above our minimum refrigerator temperature of 0.25 K, the temperature gain due to the absorption of a single photon is below 6.5 nK, which confirms our previous assumption that ∆T << T . Finally, we use the single-photon temperature gain to determine the corresponding increase in bulk resistance. Figure 5 depicts the experimental values for zero-bias resistivity as a function of temperature in bulk Cd 3 As 2 in the superconducting regime. For temperatures above 0.35 K, the resistivity steadily increases with temperature. We will therefore use 0.35 K to 0.45 K as the range of baseline temperatures for which we will determine the single-photon bulk resistance gain. For a square lattice surface, the bulk resistance scales linearly with resistivity as 1/d, where d denotes the lattice depth. Therefore, the single-photon resistance gain relates to the slope of the resistivity with respect to temperature (dρ/dT ) and the single-photon temperature gain (∆T ) as follows: For the aforementioned sample dimensions, d = 20 nm. Substituting the expression for ∆T from Eq. (48), we find that the single-photon resistance gain ∆R solely becomes a function of the baseline temperature T :  Figure 6 depicts the resistance gain due to the absorption of a single photon for baseline temperatures ranging from 0.35 K to 0.45 K for the selected photon frequencies of 5 GHz and 10 GHz. For temperatures of 0.39 K and above, the single-photon resistance gain will be greater than 1 µΩ for photon frequencies as low as 5 GHz, an increase which is certainly measurable using a commercially available micro-ohm meter (such as the Keysight 34420A NanoVolt/Micro-Ohm Meter by Keysight Technologies) or with a Corbino geometry sample which can even measure sub-micro-ohm resistance [42]. This property can therefore be exploited in order to precisely determine the number of absorbed photons for a known frequency.

VI. DISCUSSIONS AND CONCLUSION
We demonstrated a microwave photon-number resolving detector based on the topological surface states of Cd 3 As 2 material. The number of photons absorbed is produced after measuring the increased resistivity of the superconducting bulk. For this, we derived in detail how much bulk temperature would elevate as a function of the absorbed number of photons and the photon frequency. We showed that the energy of the absorbed photon is rapidly transferred firstly to the rethermalized distribution of the surface state electron band population. Then, the electron band energy is quickly transferred to the bulk phonon modes through the deformation potential coupling. The bulk temperature is thus elevated, and finally, the superconducting bulk increases resistance, which is measured to resolve the absorbed number of photons. To address how quickly the energy is transferred from the surface electron to the bulk phonon modes, we derived the deformation potential electron-phonon coupling rate by calculating the transition matrix element and the phase space volume. As a result, we concluded that the coupling time constant ranged from nanoseconds to microseconds. Therefore, it is expected that the number of absorbed photons would be measured within several milliseconds after the absorption happens.
Our proposed scheme accomplishes rapid photon detection based on quick (or even continuous) and accurate bulk resistance measurement. Direct measurement of the elevated temperature in bulk does not provide a feasible path due to the slow detection speed and the measurement noise in the extremely small differential temperature. It is essential to understand why the use of Cd 3 As 2 bulk's semimetal feature for absorbing microwave photons is avoided. Recall that, if the baseline temperature is set above the critical temperature, the bulk's electronic bands do not open a gap f ∆ as shown in equation (2), which allows the bulk electrons to be excited by the microwave photons. However, detecting the excited electron is extremely difficult for two main reasons. First, the bulk photoelectron may easily join the resistance-measuring current and be lost in the measurement process. Second, the photoelectron's energy transfer to the bulk temperature is extremely inefficient due to the reduced phase space of 3D electrons, risking the loss of photoelectrons via radiative decay rather than energy transfer to the bulk phonon modes. In contrast, the photon absorption from the surface state electrons almost surely transfers the energy to the bulk phonon modes.
Equally important is understanding the difference between our proposed scheme and an alternative device structure of a Dirac 2D material such as graphene on the surface of a bulk superconductor. A pure graphene layer indeed does not possess a superconductor state [43], and thus can be used as a Dirac cone photon absorber of microwave photons even at a very low temperature. However, it is more difficult to fabricate this device than Cd 3 As 2 which simultaneously has both bulk superconductor and surface states. In addition, the hybrid structure suffers from inefficient electronic energy transfer to the bulk phonons due to the mismatch of lattice constants. Instead, as previous research on graphene singlephoton detectors has shown, the inefficient electronic energy transfer to phonons is used for efficient capture of the photoelectron in the electrodes [44]. However, in this case, the photon-number resolving feature is lost. In comparison, our scheme utilizes the surface state electrons of Cd 3 As 2 as a microwave photon absorber and the bulk superconductor of the same material for detecting the number of photons absorbed. The distinct advantage of our method is to provide a deterministic photon-number resolving capability in microwave photon detection.
We now discuss the design strategy of maximizing the photon absorption probability of the device. Note that each crystal surface features an absorption rate of 2.3 %. Therefore, it is possible to have a near unity photon detection efficiency if one vertically stacks few hundred bulk crystal layers in a heterostructure (such that they are in series from the point of view of the incoming photon), while measuring the bulk zero-bias resistivity for each of the crystals separately. With the advent of more advanced manufacturing technique, such heterostructure is increasingly becoming possible [45].

ACKNOWLEDGMENTS
Sandia National Laboratories is a multimission laboratory managed and operated by National Technology & Engineering Solutions of Sandia, LLC, a wholly owned subsidiary of Honeywell International Inc., for the U.S. Department of Energys National Nuclear Security Administration under contract DE-NA0003525.

Appendix A: Matrix Elements for Interaction Between Surface Electrons and Bulk Phonons
In this section, we estimate the matrix elements corresponding to interaction between the surface electron modes and bulk phonon modes by building from the analogous elements for the bulk electron-phonon interaction. Labeling the direction perpendicular to the surface asẑ, we represent the surface state wavefunction as a product of the xy-plane wavefunction and a pulse-like function of z: where ψ(z) is expressed such that its amplitude-squared becomes a broadened Dirac-delta function with a width of a: Here, a denotes the approximate width of the surface state. Next, we seek to express a surface mode as a superposition of bulk modes by decomposing ψ(z) into a superposition of modes with well-defined z-wavevector k z : where |k z denotes the plane wave state with wavevector k z , taking the following form with respect to the lattice depth d when projected onto the position space: From the Heisenberg uncertainty principle, we intuitively know that the range of z-direction momentum is approximated as ∆p z ≈ h/a, leading to a wavevector range of ∆k z ≈ 2π/a. We quantitatively determine the superposition coefficients c(k z ) as follows: For low values of k z a (i.e. k z a 1), the coefficient can be estimated as a/d. Since the span of each plane wave state in reciprocal space is 2π/d, this accords with the intuition that the overall recpirocal space in the z-direction spans a length of 2π/a, subdivided into a/d wavevectors with a roughly uniform superposition coefficient for each. The intraband matrix element corresponding to the emission of a phonon of wavevector q in branch µ by a surface electron can thus be expressed in terms of the analogous matrix elements for bulk electrons: Note that this approximation is valid specifically if the maximum amplitude of the emitted phonon wavevector is much smaller than the maximum amplitude of k z , i.e. if q max << π/a, which holds true for long-wavelength acoustic phonons in the linear dispersion regime.
As a final step in the generic matrix element calculation, we can show that the carrier-phonon matrix element is exactly invariant in the initial carrier wavevector k z , since each k z corresponds to a plane-wave state with well-defined z-momentum k z . Specifically, the ma-trix element of a spatial function f (r) connecting an initial carrier state |k to a final state |k − q (where the wavevectors are three-dimensional) is simplified as follows: As a result, the terms in the summation in Eq. (A6) are equivalent. Since the total number of valid planewave wavevectors in the summation is d/a (as previously discussed), the matrix element H emit becomes invariant in the lattice depth d, as desired: Therefore, the electron-phonon matrix element M µ,q kxy,pxy corresponding to the transition from an initial electronic wavevector k xy to a final wavevector p xy via the emission of a phonon of wavevector q in branch µ is equivalent to the bulk matrix element M µ,q k,p , where k z = 0 and p z = −q z .

Appendix B: Calculating the Surface Electron Cooling Rate
We start by representing the cooling rate of the surface electrons as the following summation over initial electron wavevectors k, final electron wavevectors p, and phonon wavevectors q, with the direction of k defined as the xaxis and the projection of p on the xy-plane labeled as p xy : Given an overall lattice volume V and surface state area A, the discrete summation over k and p can be converted to integrals as follows: 3 dp xy p xy dθ p dp z .

(B2)
Also, the Dirac delta can be re-written as a function of p xy : Therefore, the Dirac delta collapses the integral over p xy to p xy = k − v s q/v F . The rate Q reduces to the following form: Next, we apply the Kronecker delta, which restricts the summation over phonon wavevectors to q = k − p. Also, as discussed in the main text, we can apply the approx-imation p ≈ p z >> k for nearly all of the phase space. This argument is further strengthened by the fact that the integrand approaches zero for low values of q for which this approximation is the weakest. Therefore, we set q = p z . Since all terms in the integrand are now functions solely of k or p z , we have collapsed the integral over p to an integral over p z multiplied by the circumference of the circle generated by making a plane cut (along the p x p y plane) through the cone. Also, the fact that the radius of the cone at a generic value of p z equals k − (v s /v F )p z implies that the circumference of the aforementioned circle is 2π(k − (v s /v F )q z ) and the upper bound of q z (corresponding to the cone vertex) is (v F /v s )k. Therefore, Q simplifies to the following: Note that the factor of 2 in front derives from the fact that the available phase space is actually a double cone, extending in both the + and − directions along the p zaxis.
We now examine the parenthetical term corresponding to the net phonon number. In the limit ∆T << T e , T L , where ∆T = T e − T L , we solve for a generic value of T ≈ T e ≈ T L to first order in ∆T : (B6) The exponential term will set an approximate upper bound for the emitted phonon energy. For low temperatures (e.g. T = 0.45 K), this ensures that the emitted phonons fall within the long-wavelength acoustic mode limit.
Next, we write out the expression for the electron occupation number difference between the initial and final states: We are now in a position to simplify the double integral by defining the variables x = v F k/(k B T ) and y = v s p z /(k B T ), yielding the following expression for Q: dy(x − y)y 3 e y (e y − 1) 2 1 e x + 1 − 1 e x−y + 1 .

(B8)
Note that our analysis so far has been restricted to the case of intraband carrier transitions through carrierphonon interaction. Unlike the case of bulk carriers, where the band structure prohibits interband carrierphonon scattering, the carrier dispersion for the surface state allows such scattering. The phase space area for this interaction can be constructed by merging the coni-cal phase space for the electron-phonon interaction with the corresponding cone for the hole-phonon interaction, with the pair of cones meeting at p z = v F k/v s . Just like we restricted the range of values for p z for the intraband scattering to p z < v F k/v s , we will restrict the range for interband scattering to p z > v F k/v s . Recall that the amplitude-squared of the carrier-phonon matrix element for Cd 3 As 2 is identical for interband and intraband interaction except for the following proportionality [25]: where θ represents the angle between k and p, and s = 1, −1 for intraband and interband interactions, respectively. As mentioned previously, p z >> k for nearly all of the phase space area for the intraband interac-tions, and since the interband interactions involve even higher values for p z , this is clearly true for such processes as well. Therefore, for both intraband and interband interactions, we can make the approximation cos θ ≈ 0, causing the interband and intraband matrix element amplitude-squared values to scale linearly with the phonon wavevector in the same manner. In order to construct the equivalent of Eq. (B8) for interband carrierphonon scattering, the only changes we make are thus the range of p z , as previously mentioned, as well as the term in the integrand corresponding to the cone radius, which flips as k − v s p z /v F → v s p z /v F − k: