Coulomb blockade thermometry beyond the universal regime

The charge localization of single electrons on mesoscopic metallic islands leads to a suppression of the electrical current, known as the Coulomb blockade. When this correction is small, it enables primary electron thermometry, as it was first demonstrated by Pekola et al. (Phys. Rev. Letters, 73, 2903 [1994]). However, in the low temperature limit, random charge offsets influence the conductance and limit the universal behavior of a single metallic island. In this work, we numerically investigate the conductance of a junction array, and demonstrate the extension of the primary regime for large arrays, even when the variations in the device parameters are taken into account. We find that our simulations agree well with measured conductance traces in the submillikelvin electron temperature regime.


I. INTRODUCTION
Single electron charging effects are prominent in isolated nanoscale conductors below a characteristic temperature scale k B T e 2 /2C Σ [1,2], where k B is the Boltzmann constant, e is the electron charge and C Σ is the total capacitance of the metallic island. In this low temperature regime, charge sensing [3] and direct current measurements through tunneling contacts with low conductance, G e 2 /h [4] yield a well-defined number of excess electrons on small metallic islands and an exponentially suppressed tunneling current. Increasing the temperature to k B T > e 2 /2C Σ results in the breakdown of Coulomb blockade, and it was shown that in this regime, the conductance suppression is universal. This enables a primary measurement of the electron temperature [5], which is insensitive to the device geometry, the external magnetic field [6] and the electrostatic environment of the device [7]. This property makes Coulomb blockade thermometry favorable for metrological applications [8][9][10][11], and several experiments were performed in cryogenic temperatures ranging down to the millikelvin [12][13][14][15][16][17][18] and submillikelvin regime [19] as well as up to 60 K [20].
The practical upper temperature range of a Coulomb blockade thermometer (CBT) is limited by the readout of the low bias conductance suppression, which is inversely proportional to the temperature [10,20]. In the low temperature limit, the conductance depends on the electrostatic environment of the device, which imposes an effective offset charge [4,21]. This offset charge depends on the random, uncontrolled population of charge traps at nearby interfaces and charged impurities inside dielectrics [22,23] leading to a statistical error for primary electron thermometry by CBTs.
Precision electron thermometry requires a characterization of this error source. In particular, the recent advent of ultralow-temperature quantum electronics [24][25][26] requires electron thermometry spanning several orders of magnitude in temperature. Here, we use a statistical Monte Carlo approach to investigate the temperature error sources of CBTs in the low temperature regime, and demonstrate that this numerical procedure provides a unified description of a CBT, matching the predictions of the commonly used single island master equation (ME) model [5] in the high temperature, universal regime. Importantly for Coulomb blockade thermometry, we show that the effect of random offset charges is suppressed in the case of large arrays, which enables reliable temperature measurements deeper in the low temperature regime, even in the presence of variations in the device parameters. Finally, we demonstrate the applicability of our numerical results by a direct comparison with ultralow-temperature experimental data, where previous numerical and analytical models fail to describe the conductance of the CBT.

II. ELECTROSTATICS OF THE CBT DEVICE
We consider the device geometry outlined in Fig. 1, which displays the typical implementation of a CBT consisting of N tunnel junctions in series enclosing N −1 islands. Each junction possesses a resistance R i and a parallel capacitance C i . Crucially for the current work, we also consider the offset-charge of each island, q o,i = −C o,i V o,i yielding a total charge q i = en i + q o,i , where the number of excess electrons on the island is n i . Throughout this work, we assume that the effective gating capacitance C o,i is much smaller than C i and C i+1 for all islands and investigate the applicability of this limit for experiments in Section 7.
The array is symmetrically biased by two voltage sources, fixing the electrostatic potential on the leftmost and rightmost nodes at −V /2 and V /2 respectively, leading to a total voltage drop of V . We neglect environment-assisted tunneling processes [27] assuming a low-impedance electromagnetic environment surrounding the device.
In our analytical and numerical calculations, we always consider a single series of junctions first and sum up the resulting conductances to model N × M arrays. This method is equivalent to neglecting the capacitive coupling between islands in neighboring chains, which is justified in typical device geometries.
We now consider the electrostatic potential ϕ i of each island for i = 1 . . . N − 1, which allows us to express the total charge on the island: where ϕ 0 = −V /2 and ϕ N = V /2. We can write this equation in a matrix form, C · ϕ = −q with where C Σi = C i + C i+1 , the total capacitance of island i. Finally, the electrostatic potentials define the voltage drop over each tunnel junction as V i = ϕ i+1 − ϕ i and the total electrostatic energy is given by We now consider a homogeneous array in the absence of offset charge, where all C i = C and q o,i = q o = 0, at zero bias voltage, V = 0. In this limit, the electrostatic energy of a single electron on island k is calculated using the net capacitance C(k −1 + (N − k) −1 ), which yields Substituting C Σ = 2C and k = 1 or k = N − 1, Eq. (4) yields the charging energy expression which defines the dimensionless inverted temperature u = 2E C /k B T [5,28].
Turning to the effect of a finite offset charge q o,i = q o on all islands of a homogeneous array, we derive the condition where the change of F (see Eq. (3)) is zero in response to a single electron tunneling event (isoenergetic tunneling). With zero excess electrons on all islands, charge conservation yields the following expression for the voltage drops over the tunnel junctions: and i V i = V . If we now insert a single excess electron onto island k, V k+1 − V k = −(q o − e)/C holds with the other differences remaining the same. We can now zero out ∆F ∝ i V 2 i − i V 2 i with the following index shift: which is fulfilled for any k when q o = e/N , yielding ∆F = 0 independent of the bias voltage V = i V i . This result has important consequences to the charge degeneracy of the device. The well-known charge degeneracy condition of q o = e/2 in case of a single island changes to q o = e/N for an array of N tunnel junctions, where the total electrostatic energy of an empty chain is equal to the energy when a single excess electron is present on any of the N − 1 islands. We illustrate this effect by plotting F at V = 0 as a function of q o in the case of 0 and ±e excess charge in the case of N = 2 ( Fig. 2a) and N = 36 (Fig. 2b), where different positions of the excess charge k = 1 . . . 4 are shown, all intersecting at q o = e/N . This result shows that in the low temperature regime, we should expect the highest conductance of a CBT array at this offset charge corresponding to minimal Coulomb blockade [29].

III. NUMERICAL MODEL
We now introduce the numerical method leading to the voltage-dependent differential conductance G(V ) of a CBT array at an arbitrary temperature. We base our calculations on the Markov chain Monte Carlo (MCMC) method [30,31], which allows us to investigate arbitrary array sizes, offset-charge configurations, and variations in the parameters of the tunnel junctions.
The first (last) N − 1 columns describe the forward (backward) tunneling of one charge quantum. Next we evaluate the resulting potentials in a concise form ϕ = C −1 Q and calculate the associated free energy change for a single electron hopping through junction i [30]: Here, ϕ 0 = ϕ 0 = −V /2 and ϕ N = ϕ N = V /2 at a bias voltage V (Fig. 1b). We note that for a homogeneous array with C k = C, the above expressions lead to ∆F i = e/2(ϕ i±1 + ϕ i±1 − ϕ i − ϕ i ), which was also used to find the analytic high temperature G(V ) curve for the N = 2 case [5]. Based on ∆F i , we calculate the forward and backward tunneling rates through junction i with a resistance R i : The tunneling current I = ±eΓ ± i is a result of the single charge tunneling events. Following standard Monte Carlo methods, we consider the current charge configuration q and evaluate the corresponding Γ ± i values. Then, we randomly select the realized tunneling process with a probability proportional to the corresponding tunneling rate. In addition, we follow the variance reduction method introduced by [31,32] to calculate I. Here, we make use of the current conservation I = I 1 = I 2 = ... = I N along the tunnel junction chain, and express the total current as I = ∆Q/∆t using [32]: where the Monte Carlo step p runs in the range of 1 . . . P , which we set to reach a relative error less than 10 −4 , typically achieved with P < 10 5 . We numerically acquire the differential conductance as G(V ) = (I(V +δV )−I(V −δV ))/2δV , where δV k B T /e. Finally, we recover the asymptotic tunnel conductance at high bias voltages as G −1 First, we discuss the dimensionless zero bias conductance G(0)/G t , which we plot in Fig. 3 for various N values and for the two limiting cases in homogeneous offset charge, q o = 0 (dashed lines) and q o = e/N (solid lines). We benchmark our numerical MCMC results in the well-studied high temperature regime, where u 1. Here, we recover the earlier results [5,29] that the conductance collapses onto a single curve independent of N and q o . The conductance in this limit has been calculated earlier by writing a master equation for forward and backward tunneling processes for a single island device with N = 2. Here, the conductance suppression ∆G = G t − G(0) follows a series expansion in u 1 [5,29]: In the low temperature limit, we find that the conductance depends on q o , and the q o = 0 case corresponding to strong Coulomb blockade results in an exponential suppression of the conductance. However, at minimal blockade, q o = e/N , our results yield a low temperature saturation conductance corresponding to G(0)/G t = 1/N (see inset of Fig. 3), which is the extension of the N = 2 case discussed earlier [29].
We understand this result based on the charging energy expression summarized in Fig. 2, where all N − 1 possible states with a single excess electron and the empty chain have the same electrostatic energy if q o = 1/N . Since all other charge configurations are higher in energy, only these N states are occupied at zero temperature, with an equal probability of P = 1/N . Taking the T → 0 limit of Eq. (10), we find if ∆F > 0 (forward tunneling) and Γ = 0 otherwise. We calculate the current as I = PΓ through any of the tunnel junctions. At low bias, where eV E C , we get ∆F = eV /N , which yields I = V /(RN 2 ). With G = dI/dV and G t = 1/(N R), we arrive to the dimensionless conductance corresponding to the numerically acquired result. This conductance and its exponentially suppressed counterpart at q o = 0 demonstrates the increasing sensitivity of a CBT to the offset charge as the temperature is lowered below E C . Furthermore, in this regime, the commonly used analytical master equation model with N = 2 fails to predict G(0)/G t regardless of q o , underlining the importance of MCMC numerics for low temperature Coulomb blockade thermometry. Next, we perform the analysis of the finite bias G(V )/G t curves, which is commonly used for primary calibration in the universal regime. Here, the full width at half minimum [5,28] with ζ ≈ 5.439 enables primary thermometry if ∆G/G t 1 or the characterization of E C for secondary thermometry with G(0)/G t using Eq. (13). We first demonstrate that MCMC calculations reproduce the master equation conductance traces in the universal regime and yield the same V 1/2 (Fig. 4).
In the low temperature regime, we observe that G(V ) depends on the offset charge q o (Fig. 5). The q o = 0 case displayed in Fig. 5a was discussed earlier in the context of charge soliton propagation on tunnel junction chains, where the threshold voltage was found to be eV t = N E C at 1/u → 0 [33,34], in agreement with our numerical results. In contrast, the q o = e/N case shown in Fig. 5b exhibits a gapless behavior, with no well-defined threshold voltage. However, in agreement with the zero bias conductance data displayed in Fig. 3, a conductance suppression is observed which reaches G(0)/G t = 1/N in the zero temperature limit.

IV. THE ROLE OF RANDOM OFFSET CHARGES
Thus far, we considered a uniform offset charge q o applied on all islands of the CBT. In the absence of externally controlled gate electrodes, the offset charge is randomized by the electrostatic environment [35], which is typically  described as an ensemble of two-level charge fluctuators embedded in tunnel barriers or material interfaces [22,23,36].
Here, we assume that any rearrangement of offset charges happens much slower than the conductance measurement timescale, which is in the order of a second, so that we can consider a fixed set {q o,i } for a MCMC calculation run, and can average the resulting conductance values to acquire the offset charge-averaged conductance at a given temperature and bias voltage.
Typical charge trap densities exceed 10 4 µm −2 on Si-SiO 2 interfaces [37] and in the oxidized surfaces of metal thin layers [36], which, together with a typical metallic surfaces areas exceeding 1 µm 2 of a single island [18], results in offset charges many orders of magnitude larger than e. In a CBT device, where the islands are tunnel coupled to electrodes, this charge will be reduced by tunneling until a stable offset charge configuration is reached. In the single island case (N = 2) this results in an offset charge uniformly distributed in the range of −e/2 . . . e/2, since any offset charge outside of this interval would enable a free energy decrease by tunneling of a charge ±e. While the same range is often assumed in the case of long arrays, the system can be trapped in a local minimum of the electrostatic energy with tunneling events towards the global energy minimum being inhibited over realistic timescales in the low temperature regime.
We demonstrate the presence of such charge configurations by initializing arrays of different N with a random q o uniformly distributed with a width of ±100e, and a stable offset charge configuration is calculated by minimizing the free energy by a Markov chain of randomly selected tunneling events. The final charge configuration is recorded when all possible tunneling events out of this state increase the electrostatic free energy of the array. We show the resulting probability densities in Fig. 6 with taking the N = 2 case as a reference. Notably, the resulting offset charges always fall into the ±e window, however we find an increasing probability of being outside of ±e/2 with increasing array size. Based on this numerical result, we follow the above initialization procedure for all MCMC calculations, instead of selecting the offset charge uniformly between ±e/2 as it was done earlier [38].
We showcase the statistical variations of G(0)/G t due to the random offset charges in Fig. 7a for different array lengths. Remarkably, the resulting distribution narrows with increasing N . This is in contrast with the uniform offset charge case shown in Fig. 3, where the difference between the maximum and minimum conductance increases as a function of N owing to the strong exponential suppression of fully blockaded arrays.
We note that the width of the conductance distribution further decreases in the typical CBT geometry, where several junction series are connected in parallel in order to increase the device conductance. Assuming an independent offset charge distribution for all arrays, the standard deviation of the conductance is expected to decrease with 1/ √ M for an N × M CBT device. We demonstrate this effect by plotting the probability density functions of the sum of M = 50 randomized conductance values in Fig. 7b.
Both for M = 1 and M 1, we also observe a saturating distribution when increasing N . This defines the long array limit, where the statistical temperature error of the CBT does not depend on the array length. To simplify further analysis and comparison with our previously acquired experimental data [18,19], we use N = 36 henceforth, which falls in the long array limit in our temperature regime of interest, 1/u > 0.1. We also note that several other experiments in the millikelvin regime were performed with CBTs of a similar or larger array lengths [14,15,17].
Next, we evaluate the offset-charge averaged conductance as a function of the temperature in a demonstration of Coulomb blockade thermometry beyond the universal temperature regime. Here, we plot the expectation value of G(0)/G t as a function of the dimensionless temperature 1/u (Fig. 8). Remarkably, the confidence intervals corresponding to ±3σ remain much narrower than the full conductance window defined by q o = 0 (dashed line) and q o = 1/N (solid line). Importantly for thermometry applications, we can evaluate the resulting temperature error ∆T , and find that even for M = 1, reliable temperature measurements with ∆T /T < 0.05 are possible if 1/u > 0.2, in contrast with the 1/u > 0.4 limit calculated earlier for N = 2 [29]. Furthermore, CBT devices with several junction series in parallel are even more suitable in the low temperature regime: a realistic CBT with M = 100 reaches the same temperature error at 1/u = 0.1 demonstrating the possibility of precision low temperature electron thermometry by using scalable fabrication of large arrays.
We display the calculated finite bias conductance curves of the same device in Fig. 9, which is the low temperature extension of the results presented in Fig. 4. Remarkably, the full conductance traces remain sensitive to the temperature even when 1/u 1, which allows for a primary calibration of the CBT in the low temperature regime, despite the absence of offset-charge control on the individual islands.

V. DISORDER IN THE TUNNEL JUNCTION PARAMETERS
Previously we argued that the effect of random offset charges can be mitigated by using devices with large N and M . However, for realistic device modeling, we also have to consider the non-uniformity in the tunnel junction resistance R i and island coupling capacitance C i values. These variations can be due to lithographic inaccuracies as well as thin film thickness variations. Previous studies addressed this error source assuming that the junction area inhomogenity of in-situ created AlO x tunnel barriers dominates [7,28], equivalent to a constant R i C i . However, recent studies demonstrated the spatial variations in the oxide layer thickness [39] and new fabrication techniques utilizing ex-situ via tunnel junctions were developed [18,40], requiring a separate evaluation of the effects of the capacitance and resistance variations on the CBT accuracy.
We first discuss the effect of the capacitance disorder by calculating a conductance histogram of an N = 36 and M = 5 device with a uniform distribution of junction capacitances, C ±∆C/2 (Fig. 10a). Here, ∆C/C = 0 corresponds to a homogeneous array, where the conductance distribution is governed by the random offset charges. Increasing ∆C/C results in a broadening of the conductance histogram as well as a quadratic drop in the expectation value of G(0)/G t , which translates to a decrease in the inferred temperature (see Fig. 8). We note that we kept 1/u constant for different ∆C values, using the expectation value C in Eq. (5). We show in Fig. 10b that the impact of ∆C/C on thermometry depends on the temperature by displaying the relative temperature error ∆T /T defined as the 3σ confidence intervals as a function of ∆C/C at different 1/u values. Our results indicate that in the low temperature regime, 1/u < 0.4 and realistic capacitance variations of a few percent, the measurement error is dominated by the randomness of the offset charges.
Next, we evaluate the case of variations in the individual tunnel junction resistances, and perform the MCMC simulations with a uniform distribution of R with a width of ∆R (Fig. 11). The resistance variations cause a weak quadratic decrease in the expectation value of the temperature, however even the extreme case of ∆R/R 1 yields a smaller systematic and statistical error than the capacitance disorder in the ∆C/C ∼ 0.1 regime. When comparing the N = 2 and M = 1 (Fig. 11a) single island CBT with the N = 36 and M = 5 case (Fig. 11b), we find a marked reduction in the statistical temperature error for the long array, further attesting to the importance of large array devices for low-temperature thermometry.  9. The offset charge-averaged normalized differential conductance curves of a N = 36 device as the function of the dimensionless bias voltage in the low temperature regime where u > 1 calculated by the MCMC model. Note that the voltage scaling is the same as in Fig. 4.

VI. COMPARISON WITH EXPERIMENTAL DATA
We now demonstrate the correspondence of the MCMC simulations with experimental data taken in the submillikelvin electron temperature regime. We fabricated our Coulomb blockade thermometers with ex-situ Al/AlO x /Al via tunnel junctions connecting adjacent Al islands separated by sputtered SiO x interlayer dielectric [18,40]. This process, in contrast to the frequently used Dolan bridge technique [41], enables an independent tuning of the tunnel junction area and the capacitive coupling between the islands. Crucially for the applicability of our numerical model, our device fulfills the C o C requirement with a large fraction of the island areas overlapping (see the inset of Fig. 12). Our device features N = 36 tunnel junctions in series and M = 5 rows in parallel. Electron cooling was performed by on-and off-chip nuclear refrigeration [16], using large volume electrodeposited indium fins. When the nuclear demagnetization was carried out on the device, we observed [19] that the measured electron temperature scales with the effective magnetic field magnetic field corresponding to the quadrupolar crystal field in indium, consistent with prior literature data [42]. We found that the distortion of G(V ) because of self-heating [43] was negligible. The details of the device fabrication and measurement setup are published elsewhere [18,19].
First, we measured our CBT in the universal regime, where the ME model provides a reliable calibration of the island capacitance, C = 670 ± 2 fF, which we will use as a fixed parameter for the rest of the analysis. In the universal regime, where G(0)/G t > 0.8, we can get reliable fits using the ME model (see the two uppermost curves in Fig. 12). Below this threshold, however, we first fit the temperature comparing the zero bias conductance with the MCMC model prediction (see Fig. 8), and then confirm the validity of the fit in the full voltage bias range (four lower curves in Fig. 12). With this procedure, we get an excellent agreement between our numerical simulations and experimental data down to T = 470 µK corresponding to 1/u = 0.175, where we estimate a relative temperature error of ∆T /T = 0.049 (see Fig. 8). This result demonstrates that accurate Coulomb blockade thermometry is possible outside of the commonly considered universal regime.
Finally, we discuss the temperature error associated with the variations in the tunnel junction parameters and compare the resulting temperature error with the previous result. Based on the maximum alignment error of 200 nm of the laser lithography used to fabricate our device, we infer a relative error in overlap and in the island capacitance ∆C/C = 0.024 (see vertical arrow in Fig. 10b), which in our temperature of interest contributes with less than 0.01 to the relative temperature error. We characterized our tunnel junctions to exhibit an areal resistance 12.8 ± 0.8 kΩµm 2 [18], corresponding to ∆R/R = 0.063 (see vertical arrow in Fig. 11b), which has a negligible effect on the temperature error. Our analysis attests to the insensitivity of CBTs to realistic fabrication inaccuracies in the low temperature regime, in agreement with similar numerical studies in the universal regime [5,7,20,32].

VII. CORRECTIONS WITH A FINITE GATING CAPACITANCE
We now turn to the case where the island gating capacitance C o,i is finite, but still much less than the C i capacitances between the islands (Fig. 1). First, we establish the validity of this limit by estimating the self-capacitance of the islands in our experiments [18,19], which is assumed to be dominated by the indium cooling fins [18] with linear dimensions of 25 × 50 × 140 µm 3 . To provide an estimate, we use the self-capacitance C o = 4πε 0 d · 0.66 of a cube [44] with d = 50 µm side yielding C o ≈ 3.7 fF and C o /C ≈ 5.5 × 10 −3 .
We use the numerical model discussed in Section 3 using the capacitance matrix C (see Eq. [2]) with a modified C Σi = C i + C i+1 + C o,i . As we demonstrated in Fig. 10  caused by the random offset charges, therefore we evaluate the offset-charge averaged zero bias conductance curves as a function of temperature for finite C o /C ratios while keeping C Σ constant with no disorder in the capacitance values and tunnel junction resistances. Similarly to our other calculations, we set N = 36 to enable a direct comparison with our experiments.
We summarize our calculations in Fig. 13, where we display the temperature dependence of the normalized zero bias conductance (Fig. 13a) for a C o /C ratio of 5.5 × 10 −3 and 2 × 10 −2 , which we compare to the C o /C = 0 model also shown in Fig. 8. We recover the high temperature, universal regime, where the curves merge, however we find an increasing systematic deviation when lowering the temperature to the the intermediate [29] and to the low temperature limit. In Fig. 13b, we characterize the resulting relative temperature error resulting from using the C o = 0 model in case of these two finite C o /C ratios. We note that ∆T /T > 0, thus neglecting C o yields an overestimate of the device temperature. Our numerical results show that this error is of the order of 1% in our temperature regime of interest, and it increases with the C o /C ratio of the device, demonstrating the importance of the specific CBT device design for ulta-low temperature thermometry.
We conclude this analysis by comparing this systematic deviation to the statistical uncertainty caused by the random offset charges for M = 5, relevant for our device. We find that in the low temperature limit 1/u < 0.2, the statistical uncertainty dominates. This is in contrast to the intermediate regime 1/u 1, where the finite C o /C ratio plays a more important role. Our numerical results demonstrate that CBT devices with large area overlaps between neighboring islands [15,18] are preferred in the intermediate regime to devices with shadow-evaporated tunnel junctions [20] and large surface area cooling fins [16] featuring higher C o /C ratios.

VIII. CONCLUSIONS
Using Markov-chain Monte Carlo simulations, we demonstrated that Coulomb blockade thermometry can be extended beyond the universal regime. Our results show that long arrays with many islands in series benefit from the narrowing total conductance distribution even in the case of random offset charges, and the statistical temperature errors are further reduced by adding several junctions series in parallel, which is a typical geometry for practical CBTs. These contributions lead up to a factor of four increase in the usable temperature range of a CBT with a realistic number of tunnel junctions, impacting precision thermometry applications. We demonstrated the applicability of the numerical results to experimental data taken in the sub-millikelvin temperature regime. These temperatures fall outside of the universal regime for a CBT featuring a charging energy exceeding 1 mK, which is necessary to perform Comparison of ultralow-temperature experimental data and theoretical curves using a device with N = 36. The island capacitance of C = 670 ± 2 fF is obtained from fitting against the master equation (ME) model in the universal regime, see the uppermost two curves. The MCMC fit curves at lower four temperature values are acquired by fitting the conductance at zero bias, and then simulating the rest of the curve with no additional parameters. In this regime, the ME model does not describe the conductance curves. For all sets, solid lines display the calculated curves and dots represent the experimental data taken at a fixed magnetic field of 50 mT. The inset shows the schematic cross-section of the measured CBT featuring via tunnel junctions between neighboring islands (dark grey) through the dielectric layer (light grey) yielding C = CJ + Ce, where the capacitance contribution Ce from the overlapping islands dominates.
a primary calibration at electron temperatures in the range of 10 . . . 100 mK attainable in commercial dilution refrigerators. Finally, we showed that at sub-millikelvin temperatures, device designs yielding a low gate capacitance ratio C o /C 10 −3 are necessary to suppress systematic temperature errors stemming from a finite C o .

DATA AVAILABILITY
Raw datasets and computer code are available at the Zenodo repository [45].