### Instrument design

The instrument comprises four differential vacuum stages with decreasing pressure (Fig. 1a, left). The first stage houses an electrospray ionization (ESI) source at atmospheric pressure that generates mostly desolvated charged bacteria cells. The bacteria are then immediately attracted by the second stage, a capillary at 10 mbar that is heated at 150 °C to fully desolvate the microdroplets and to prevent sticking of bacteria to the internal wall of the capillary. The charged particles at the exit of the heated capillary come into the aerolens system comprising a long capillary followed by two consecutive skimmers^{37}. Numerical simulations by the finite element method (FEM) show that the cells slow-down in the capillary from supersonic speeds (600–700 m/s) up to tens of m/s, achieving laminar flow (Fig. 1a, right). The skimmers give rise to a collimated narrow cell beam in the fourth differential vacuum chamber at 0.1 mbar that houses the resonator chip just 14 mm below. Near the resonator, bacterial cells are increasingly deaccelerated up to they ‘softly’ land on the resonator surface^{38}. Soft-landing results into well-preserved cell morphologies with small contact areas with the resonator surface (Fig. 1b). In addition, the vacuum pressure in the detector chamber is more than four orders of magnitude higher than in any class of native mass spectrometry, low enough for removing the water content of the bacteria, while high enough for maintaining the native structure of the cells.

We characterized the detection efficiency of our instrument, which can be expressed as the product of the transmission efficiency and the capture efficiency. The first quantity is the fraction of the number of particles nebulized from the solvent to ambient conditions that reaches the detector vacuum stage. The transmission efficiency is independent of the mechanical detector, and it solely depends on the instrument architecture. The capture efficiency, in turn, is just determined by the ratio of the effective sensing area of the resonator to the cross-section of the charged cell beam at the resonator’s position. The cell distribution at the resonator position after nebulization was characterized by dark-field microscopy and image processing algorithms (Fig. 1c). The cells were located in a region with full width at half maximum (FWHM) of ~1.6 × 1.2 mm, as shown in Fig. 1c. for *S. epidermidis*. The data for *E. coli* follows a very similar behavior (Supplementary Fig. 1). We estimate a transmission efficiency of about 0.03% that is similar to previous state-of-the-art^{29,30}. Most of the lost cells arise from the distance between the ESI needle and the heated capillary. The distance must be close enough to ensure high transmission, but also sufficiently large to avoid that water droplets enter in the vacuum system. The capture efficiency as a function of the detector area is estimated by integrating the cell surface density over a rectangular region that maximizes the number of collected cells (Fig. 1d, top graph). The capture efficiency is of about 5% for the ultrathin silicon nitride membranes used in our work, which is two orders of magnitude higher than in previous NMS works operating in vacuum^{29,30}. We emphasize that the extraordinary large capture area does not compromise the mass resolution for characterizing bacteria due to the nanoscale thickness of the membranes and the narrow linewidth of the membrane vibration modes as shown below. We calculate the expected mean time between cell adsorption events as a function of the resonator area (Fig. 1d, bottom graph). In our experimental conditions (10^{9} cells/mL and flow rate of 300 nL/min), we predict a mean interval time between events of 2–3 s, which is similar to that experimentally observed. Beam-like resonators traditionally used for mass sensing, rarely exceeds capture areas of about 1000 µm^{2}, giving mean interval time between events of minutes, and even hours for the smallest devices.

### Device design and operation

We designed ultrathin 50 nm thick silicon nitride membranes with side lengths *L*_{x} = 400 *μm* and *L*_{y} = 350 μm (Fig. 2a). The small difference between the lengths, of 12.5%, plays a key role in the performance of our technology as discussed below. The vibration mode shapes of the membranes were experimentally measured by stroboscopic digital holographic microscopy (Fig. 2b)^{39}. The experimental vibration mode shapes are well-described by linear elasticity theory of plates for the regime where the tension-induced stiffness of the membrane is much larger than the bending stiffness^{40}. In this case, the eigenmodes are described by \({\psi }_{{ij}}=2{\sin }\left[\frac{i\pi }{{L}_{x}}x\right]{\sin }\left[\frac{j\pi }{{L}_{y}}y\right]\), where the mode identifier (*i*,*j*) refers to the number of antinodes in the *x* and *y* directions along the long and short sides of the membrane, respectively. The membrane displacement was measured by the laser-beam deflection method (Fig. 1a and Supplementary Fig. 2)^{41}. In order to maximize the mass resolution by the inverse problem method, we simultaneously measured the frequency of six vibration modes: (1,1), (2,1), (1,2), (3,1), (1,3), and (3,2), by digital phase-locked loops (PLLs). The membrane oscillation was driven by a piezoelectric actuator just below the membrane chip. The fourth vibration mode, (2,2), was excluded because it could not be efficiently excited in our experimental conditions. The frequency spectra of the selected vibration modes are shown in Fig. 2c. The resonance frequencies are located at 470, 715, 770, 995, 1105 and 1170. These values approximately match the theoretical values given by \({\nu }_{{ij}}=\frac{1}{2}\sqrt{\left({\left(\frac{i}{{L}_{x}}\right)}^{2}+{\left(\frac{j}{{L}_{y}}\right)}^{2}\right)\frac{\sigma }{\rho }}\), for a tensile stress *σ* = 0.19 ± 0.02 GPa, assuming a density *ρ* = 2900 Kg m^{−3} (see Methods section). The corresponding Q-factors are very high in our low-vacuum conditions, of about 1500, 2700, 3000, 3000, 3500, and 4100, respectively. The frequency noise of the six vibration modes was estimated by calculating the Allan deviation^{42,43} (Fig. 2d). At short integration times, the Allan deviation decreases with the root square of the integration time and it is inversely proportional to the amplitude, which is consistent with frequency stability dominated by white phase noise. For integration times between 0.3 and 2 s, the Allan deviations reach the minimum of ∼10^{−7}, and then they increase up to achieve an asymptotic regime where the Allan deviations linearly increase with the integration time due to drift-like sources of noise. The main source of drift arises from the effect of temperature on the stress of the membrane^{44}. We minimize the eigenfrequency drift by waiting 3–4 h for thermal stabilization of the system before the measurements. The Allan deviation analysis provides a minimal detectable mass of ~0.7 fg for a particle located at the center of the membrane when measuring the fundamental resonance frequency.

We now turn our attention to the effect of the rectangularity of the membrane on the technique performance. In perfectly squared membranes, the vibration modes *i* ≠ *j*, are degenerate, i.e., modes (*i*,*j*) and (*j*,*i*) vibrate at the same frequency^{45}. However, unavoidable small imperfections due to the fabrication process or small adsorbates on the membrane break the degeneration of these modes, splitting the degenerate resonance frequency in two very close resonance frequencies^{46}. This configuration is highly unstable as any adsorption event can modify the frequency order of these two modes, leading to erroneous solutions by the inverse problem method. The rectangularity used here is high enough for stable breaking of the frequency degeneration, i.e., the vibration mode (*i*,*j*) has a lower frequency than the (*j*,*i*) for *i* > *j*, irrespective of the accumulated adsorption events. Very importantly, the rectangularity is sufficiently small to keep a quasi-symmetric capture area that is needed to optimize the detection efficiency. Moreover, the frequency density of vibration modes is significantly increased as each degenerate mode splits into two close, but uncoupled modes. Note that the frequency ratio between the 7th vibration mode and the fundamental mode is only of about 2.5, which significantly simplifies the frequency tracking of a high number of vibration modes as performed here.

### Instrument operation and data analysis

We here describe the pipeline for determination of the dry mass of bacterial cells. The process starts by measuring the fractional frequency changes of the selected vibration modes of the membrane during nebulization of the bacteria solution. Figure 3a shows a short real-time record of these signals during nebulization of *S. epidermidis* cells. These curves are characterized by plateaus interspersed with quasi-instantaneous variations, referred here to as ‘frequency jumps’, that simultaneously occur in several vibration modes. The frequency jumps are identified as “tentative” cell adsorption events if the frequency change of at least three vibration modes is at least twice the standard deviation of their corresponding frequency fluctuations. In this experiment, 11 tentative adsorption events were identified in the 27 s time span. We chose a PLL integration time of 140 ms for measuring the resonance frequencies of the membrane, which provides a good compromise between accuracy and capability for resolving temporally close frequency jumps (see Fig. 2d). For highly tensioned membranes, particle adsorption induces changes in the kinetic energy due the added mass of the particle, whereas the potential energy is negligibly affected by the stiffness of the particle^{34,35}. In these conditions, the relative frequency shift of a vibration mode induced by the adsorption of a particle at the membrane position (*x*_{0},*y*_{0}) is given by,

$$\frac{\triangle {\nu }_{{ij}}}{{\nu }_{{ij}}}=-\frac{m}{2M}{{\psi }_{{ij}}\left({x}_{0},{y}_{0}\right)}^{2}$$

(1)

where *M* is the mass of the membrane and *m* is the mass of the particle. Notice that due to the inherent symmetry of the membrane, there exist four symmetric adsorption positions where the changes in the eigenfrequencies are identical. For the sake of simplicity, hereinafter the analysis is restricted to a quarter of the membrane. We briefly summarize here the inverse problem method, leaving the details for Supplementary Note 1. Resolution of the mass and position of the particle requires at least the measurement of three vibration modes. In the absence of noise, the inverse problem is *well posed*, i.e., the inverse solution (*m*,*x*_{0},*y*_{0}) maps to a single forward solution \(\left(\frac{\triangle {\nu }_{11}}{{\nu }_{11}},\frac{\triangle {\nu }_{21}}{{\nu }_{21}},\frac{\triangle {\nu }_{12}}{{\nu }_{12}}\right)\). Notice that this one-to-one bidirectional relation cannot be obtained with other choice of modes. The frequency noise makes the problem ill posed, the problem is not invertible for all possible frequency outputs, and the inverse solution is retrieved in terms of probability rather than a deterministic solution. The accuracy of the inverse problem method increases with the number of measured vibration modes as done here. Each frequency jump is characterized by a mean value and a standard deviation. The probability density function (PDF) of the fractional eigenfrequency shift vector \(\hat{\Omega }=\left[\frac{\triangle {\nu }_{11}}{{\nu }_{11}},\frac{\triangle {\nu }_{21}}{{\nu }_{21}},\frac{\triangle {\nu }_{12}}{{\nu }_{12}},\frac{\triangle {\nu }_{31}}{{\nu }_{31}},\frac{\triangle {\nu }_{13}}{{\nu }_{13}},\frac{\triangle {\nu }_{32}}{{\nu }_{32}}\right]\), due to a particle adsorption is given by,

$${PDF}\left(\hat{\Omega }\right)=\frac{{{{{{{\rm{e}}}}}}}^{-\frac{\left(\hat{\Omega }-\hat{\triangle }\right){\Sigma }^{-1}{\left(\hat{\Omega }-\hat{\triangle }\right)}^{{{{{{\rm{T}}}}}}}}{2}}}{{\left(2{{{{{\rm{\pi }}}}}}\right)}^{\frac{{{{{{\rm{N}}}}}}}{2}}\sqrt{|\Sigma |}}$$

(2)

where \(\hat{\triangle }\) is the vector of the mean frequency shift, Σ is the covariance matrix and *N* is the number of vibration modes. By substituting the fractional eigenfrequency vector, \(\hat{\Omega }\), by its theoretical dependency explicited in Eq. (1), we build the joint probability density of the mass particle and adsorption position. Instead of selecting a recovery method such as the minimum-mean-squared error (MMSE), maximum likelihood (ML) or maximum a posteriori (MAP) probability^{47}, we here separately provide the probability density functions for the mass and position of the bacterial cells.

The proposed inverse problem method is applied to the frequency jumps observed in Fig. 3a. Inverse solutions are found for all the events, except for the fifth and sixth that cannot be clearly resolved due to the short temporal separation, of about 400 ms. The adsorption positions with highest probability for the remaining events are compared to the positions of the scattering spots at the membrane obtained by dark-field microscopy after the experiments (Fig. 3b). We observe that the predicted positions are very close to darkfield hotspots. The corresponding mass probability distributions peak at masses between 200 and 350 fg (Fig. 3c). These peaks significantly broaden as the particles are near the membrane edges, where the amplitudes of the eigenmodes are smaller. This is notably the case of particle 2, whose derived position significantly deviates from its real position and its broad mass distribution prevents accurate quantification of its mass. Later, it is shown how poorly resolved events are excluded of the analysis based on mathematical criteria. We compare the masses obtained by the inverse method with the sizes obtained by SEM images of the identified bacteria on the membrane (insets of Fig. 3c). We find that *S. epidermidis* cells are present at different growth stages, consistent with steady-state population of asynchronously growing cells. For instance, cell 3 and cell 9 were under division. We find a good correlation between mean cell mass estimated by the inverse problem method and the sizes estimated from the SEM images. The slope between the predicted mass and the volume estimated from the SEM images is of about 800 ± 200 Kg/m^{3}, similar to previously reported density values (Fig. 3d)^{29,48}.

### Determination of the mass heterogeneity of bacteria populations

A 42 min real-time record of the fractional frequency changes of the selected vibration modes of the membrane during ESI of *S. epidermidis* is shown in Fig. 4a. About 900 tentative adsorption events were found. Not all events can be solved by the inverse problem method. A small portion are artifactual and may arise from spurious electronics’ noise, external mechanical disturbances or errors in the jump selection method. Other tentative adsorption events are true, but the frequency jumps are not precisely measured because the low signal-to noise ratio. This often occurs when the cells adsorb near the membrane edges, as shown in the previous section. We impose an admissibility condition to the tentative adsorption events based on two selection rules: i) the highest probability position of the particle obtained by the inverse problem method must be separated from the membrane edges more than 6% of the length and ii) the Pearson correlation coefficient between the eigenmode square amplitudes at the highest probability particle position and the experimental eigenfrequency downshifts must be higher than 0.9 (Supplementary Note 1). Assuming that the model that relates the particle parameters to the vibrational properties of the membrane is correct, a true adsorption event gives rise to a correlation coefficient of 1 for a noise-free measurement and it decreases with the frequency noise. Moreover, false events will result into anomalous low linear correlations and eventually in adsorption positions near the membrane edges. After imposing these selection rules, the number of valid events is reduced to 70% of the tentatively identified events. The resulting data shows that the single cell biomass is highly variable, leading to rapidly changing mass distributions with the number of measured cells for *N* ≲ 100 (Fig. 4b, c). The statistical parameters that describe the mass distribution start to reach the asymptotic regime for *N* ≳ 300. The final dry mass distribution (*N* = 628) exhibits a broad distribution with mean and standard deviation of 153 ± 69 fg, discretized with narrow peaks suggesting mass values with higher probability. This effect may be attributed to the high mass resolution, with a median of 3 fg, which enabled to resolve variations of 2% in the mass of single *S. epidermidis* cells. Indeed, the worst the mass resolution, the broader the probability distribution associated to a single cell, resulting into smoother mass distributions, such as those observed by optical and electric techniques^{49,50}. Our mass accuracy is enough to describe the cell-to-cell variability, as shown below.

We next investigated the nature of the dry mass heterogeneity of the bacteria cells. Cell growth and cell division are both subject to stochastic processes due to noise in the molecular circuit dynamics, partition noise and environmental fluctuations^{1,2,49,50,51}. Cell size heterogeneity has been often quantified by statistical parameters, namely the coefficient of variation (CV). First attempts to relate the size distribution of steady-state populations of asynchronously growing cells to the source of stochasticity date back to the sixties in the works pioneered by Koch & Schaechter^{52,53}. Unfortunately, these works were latter ignored, probably because the techniques at that moment, mainly the Coulter method, lacked of the needed accuracy. We here revisit these theoretical concepts for interpreting our data. Briefly, the steady-state mass distribution of exponentially growing cells that are free of stochasticity (deterministic or canonical distribution) is given by \(\theta \left(m\right)=\frac{{m}_{d}}{{m}^{2}}\), when \(\frac{{m}_{d}}{2}\le m\le {m}_{d}\) and zero otherwise, where *m*_{d} is the mass at division. Stochasticity in the cell growth and division makes that the mass distribution deviates from the deterministic distribution, and it is this deviation the key property for describing the intrinsic stochasticity of the cell population. Koch & Schaechter equation describes the effect of cell division stochasticity on the mass distribution (Fig. 5a, left),

$$\theta \left(m\right)=\frac{A}{{m}^{2}}{\int }_{m}^{2m}g\left(x\right){dx}$$

(3)

where *g*(*x*) is the mass probability density of the cells at the instant of division and *A* is a normalization constant. The beauty of this approximation is that the coefficient of variation (CV) for the deterministic growth is 0.2017. Therefore, intrinsic stochasticity is reflected into an excess of CV with respect to this value (Fig. 5b). We also examine the effect of mass resolution on the mass distribution (Fig. 5a, right) and on the CV (Fig. 5b). This effect has been overlooked in the past, but may lead to erroneous interpretation of the data. We find that mass resolutions above 5% can significantly distort the mass distributions and induce a false increase of the CV excess. This effect is not significant in our technique, but it may be important in other techniques with lower size resolution.

We apply the Koch & Schaechter model assuming a normal distribution of the division mass to the steady-state distributions of *S. epidermidis* (*N* = 628) and *E. coli K-12* (*N* = 685; Fig. 5c, d). Mean dry mass of *S. epidermidis* and *E. coli* were 153 and 469 fg, respectively. The coefficients of variation were 0.46 and 0.37, respectively. Although cell size distributions are usually fitted by normal and lognormal distributions, these statistical models are not substantiated in any growth model^{50,51,53}. We find that Koch & Schaechter model fits with higher R-squared, and enables to infer information on the growth stochasticity parameters (Supplementary Fig. 3). The model provides a mean mass at the instant of division of 248 fg for *S. epidermidis* and 763 fg for *E. coli*. The CV of the division mass was 0.31 and 0.27, respectively. We remark again the relevance of measuring the dry mass with high accuracy to determine the sources of intrinsic stochasticity. Most studies on bacterial cell growth rely on real-time optical imaging of rod-like bacteria with lengths that range from 2 to 10 µm and widths of 0.2–1 µm. Cell volume is usually quantified by length assuming constant width due to the limited spatial resolution. Reported division length CV for rod-like bacteria lies between 10 and 15%^{50,54}, from 2 to 3 times smaller than the CV of the division dry mass of *E. coli* reported here, which question the use of length as a reliable size indicator. Moreover, the limitation in spatial resolution of the used optical techniques prevents the study of cell growth of spherical bacteria smaller than 1 µm, such as *S. epidermidis* cells studied here.