Intelligent non-colorimetric indicators for the perishable supply chain by non-wovens with photo-programmed thermal response

Spoiled perishable products, such as food and drugs exposed to inappropriate temperature, cause million illnesses every year. Risks range from intoxication due to pathogen-contaminated edibles, to suboptimal potency of temperature-sensitive vaccines. High-performance and low-cost indicators are needed, based on conformable materials whose properties change continuously and irreversibly depending on the experienced time-temperature profile. However, these systems can be limited by unclear reading, especially for colour-blind people, and are often difficult to be encoded with a tailored response to detect excess temperature over varying temporal profiles. Here we report on optically-programmed, non-colorimetric indicators based on nano-textured non-wovens encoded by their cross-linking degree. This combination allows a desired time-temperature response to be achieved, to address different perishable products. The devices operate by visual contrast with ambient light, which is explained by backscattering calculations for the complex fibrous material. Optical nanomaterials with photo-encoded thermal properties might establish new design rules for intelligent labels.

| SEM micrographs of SU-8 fibers spun upon dilution in acetone at different relative volumes (acetone/SU-8 volume): a, 0%, b, 3%, c, 7%, d, 9%, e, 12% and f, 16%. Additional electrospinning parameters: flow rate= 1 mL/h, needle-to-collector distance= 30 cm, electric field= 0.7 kV/cm. Scale bars, 2 µm (a-f) and 1 µm (insets a-e).    Blue line: as-spun fibers. Red line: fibers at the second step of photo-programming by an exposure dose of 305 mJ/cm 2 . The dashed lines indicates the peak at 911 cm -1 , assigned to epoxy ring vibration, and the reference peak at 1607 cm -1 , used to calculate the cross-linking degree, . Figure S6 | a, DSC thermal analysis of unprogrammed (pristine) nano-textured SU-8 non-wovens (bottom curve), and (from bottom to top) for non-wovens exposed by 0.3, 23, 43, 62, and 80 mJ/cm 2 (first photo-programming step). Regions of DSC profiles with inflection points (dots) indicative of Tg gradually move to higher temperature (dashed line, guide for the eye) upon increasing the exposure dose. Exothermic peaks are visible at temperature  50°C, possibly partially hindering the glass-transition behaviour for high exposure dose. The exothermic peaks are associated with polymerization induced by a very low amount of active initiator [1] and with thermal curing of the organic system [2]. b, DSC magnification of the glass transition temperature region for the unprogrammed sample. Tg0 13°C (highlighted by the arrow).
Figure S7 | Photographs of SU-8 non-wovens, deposited onto glass substrates covered with indium tin oxide and kept at different temperatures for 60 minutes. First photo-programming dose = 80 mJ/cm 2 . Tamb: ambient temperature, about 20°C. The logo "NEST" behind the samples becomes visible once upon temperature increasing and closely approaching Tw=38°C. Figure S8 | Quantitative analysis of the thermal gradient between sample regions unexposed (w/o-UV) and exposed (w-UV) to UV light during the second step of photo-programming (as in Figure 2a-d). First photoprogramming dose = 80 mJ/cm 2 . Second, spatially-selective photo-programming dose  250 mJ/cm 2 . Device operating at 55 °C. Inset: magnification view (interval 1-60 s).

Figure S9 |
Backscattering intensity maps at RGB reference wavelengths (450 nm, 530 nm, 600 nm). The field intensity maps are calculated by considering a plane wave illumination with unitary amplitude impinging perpendicularly to the structure, as in Figure 3. The maps are obtained averaging over polarization states and are calculated at about 1 m distance from the first layer of the structure. The corresponding albedo value for the three different wavelengths is A450=0.69, A530=0.74, A600=0.77. Figure S10. Albedo calculated by T-matrix formalism (circles) at  = 510 nm, as a function of the non-woven thickness. The albedo is calculated for a system thickness = 700 nm, 1.4 m, 2.8 m, 4.2 m, 5.6 m, and 7.0 m (namely up to ten layers). For non-absorbing structures, the albedo approaches the unity asymptotic value for thickness larger than the transport length (l*=3.4 m, dashed red lines) calculated from the similarity relation (see details below). Blue line: a simple model of the albedo based on the Lambert-Beer extinction decay, with l*=3.4 m and the albedo for the first layer is fixed at the value calculated by T-matrix (see details below). The inset is a close-up at small thickness, that shows the good agreement of the simple model to the albedo values calculated with T-matrix.
Figure S11 | a, Schematic representation of the assembly of a device consisting of a PET-encapsulated fibrous non-woven on paper. b-e. Representative photographs of the device at each fabrication step: b, black paper substrate; c, fiber deposition; d, UV photo-programming (step 1, dose: 72 mJ/cm 2 and step 2, dose: 305 mJ/cm 2 ); e. encapsulation with PET foils. f,g. Photographs of the device after thermal exposure at T= 55°C under static and bending conditions. Scale bar, 5 mm.

Light scattering in the T-matrix formalism
Light scattering calculations are carried out by considering a monochromatic incident polarized plane wave described by the complex amplitude: with unit polarization vector ̂and propagation vector =̂, where k=nkv and kv=/c. Here, n is the refractive index of the external medium that is fixed to n=1 since the fibers are in air. The total field outside the particle is the sum of the incident and scattered field. The scattered field is calculated from the equations obtained by applying the customary boundary conditions across the particle surface linking the internal and external fields. The overall scattering problem can be solved by the transition matrix (Tmatrix) method [3][4][5][6][7][8], based on the definition of a linear operator relating the incident field to the scattered field [3].
The starting point of the method is the field expansion in terms of vector spherical harmonics (VSH), i.e., the vector solutions of the Maxwell equations in a homogeneous medium that are simultaneous eigenfunctions of the angular momentum and the parity operators. The multipole expansion of the plane wave electric field and of the corresponding scattered wave is [4]: where ( ) ( , )denotes vector multipole fields that are regular at the origin of the coordinate system fixed inside the particle (Bessel multipoles), superscript p=1,2 refers to the values of the parity index, distinguishing magnetic multipoles (p=1) from electric ones (p=2), the ( ) denotes multipole fields which satisfy the radiation condition at infinity (Hankel multipoles) and they are identical to the the T-matrix under rotation of the scattering particle [8]. This is applicable to the calculation of the orientational average of all the quantities of interest for any choice of the orientational distribution function of the particles.

Aggregates of spheres.
We modeled the backscattering properties of nanofiber non-wovens with specific thickness considering each fiber composed of several spherical particles and exploiting the T-matrix of such aggregate of spheres. In fact, the T-matrix approach solves the scattering problem by an aggregate (cluster) of spheres without resorting to any approximation [9][10][11][12].
In particular, we consider a cluster of spheres numbered by the index , of radius ρ , and refractive index , whose centers lie at . The spheres may have, in principle, different refractive index and different radii. The field scattered by the whole aggregate is written as the superposition of the fields scattered by the single spheres: where = − . The field within each sphere is taken in the form: , so that it is regular everywhere inside the sphere. The amplitudes ( ) and ( ) are determined by applying the boundary conditions tot he fields across the surface of each sphere. Thus, the scattered field is given by a linear combination of multipole fields with different origins, whereas the incident field is given by a combination of multipole fields centered at the origin of the coordinates. In order to write the whole field in terms of multipole fields centered at we use the addition theorem [4]. In this way the scattered field at the surface of the -th sphere turns out to be [4,12]: where the operator H is a transfer matrix. Analogously, the incident field at the surface of the -th sphere is [4,12]: where 0 = 0 is the vector position of the origin. Applying the boundary conditions, we get for each  four equations among which the elimination of the amplitudes of the internal field can be applied. In this way we get a system of linear inhomogeneous equations: where we define: The quantities ( ) are the Mie coefficients for the -th sphere. The quantities which result from the addition theorem, describe the multiple scattering processes occurring among the spheres in the aggregate. The occurrence of these processes is indicative of inter-spheres coupling. The elements of the transfer matrix couple multipole fields both of the same and of different parity, with origin on different spheres. The formal solution to system (S9): relates the multipole amplitudes of the incident field to those of the field scattered by the whole object.
In order to define the T-matrix for the whole aggregate it is necessary to express the scattered field in terms of multipole fields with the same origin. Thanks to the addition theorem, the scattered field can be written as: which is valid in the region external to the smallest sphere including the whole aggregate. The previous equation shows that the field scattered by the whole cluster can be expanded as a series of vector multipole fields with a single origin (monocentered expansion) provided that the amplitudes are: Substituting for the amplitudes their expression (S11), we can define the T-matrix of the aggregate as: These are the key quantities needed to define the normalized scattering amplitude matrix of the aggregate [4] that encompasses all the observable features of the scattering process, as it is related to the flux of the electromagnetic energy that the particle scatters within the unit solid angle around the direction of observation ̂.
In order to calculate the T-matrix of a cluster we have to solve the system (S9) which has, in principle, infinite order. The system must be truncated to some finite order by including into the multipole expansions terms up to order lM, chosen to ensure the convergency of the calculations. For a cluster of Nc spheres this implies solving a system of order dM=2 Nc lM (lM +2), which may grow rather large. In our fiber calculations we used lM=4 that for the largest cluster, with 360 subunits, yields = 17280. Indeed, the inversion of the matrix M is responsible for most of the CPU time required for the calculation, which scales as 3 . We also note that, for our fiber sample, Nc=N1nl where N1~36 is the average number of subunits composing one layer, and nl is the number of layers. As a consequence, the computational demand increases with the cube of the number of the spheres, Nc 3 , and hence with the cube of the number of layers, nl 3 . For the largest cluster the computing request was about 12 GB of RAM per single simulation, taking about 8 hours on an INTEL2650v3 CPU-based computing cluster.

Scattering mean free path and transport length.
T-matrix calculations can be used to investigate the light transport properties in nanofiber samples, for which the single-scattering approximation is no more valid. Indeed this analytical approach, as clarified by Mishchenko [13], allows the radiative transfer equation to be derived from first principles. The need to consider the multiple-scattering processes has been stressed by Ishimaru and Kuga [14] who found that, as density increases, measurements might significantly depart from those expected for low-density dispersions. The multiple scattering processes within a dispersion can be included through the Foldy-Twersky equation for the propagation of the coherent field [15,16], excluding processes that imply cyclic paths (Twersky approximation). When the primary incident field is a plane wave, the number density is a constant, and the particles are in each other's wave zone an interesting solution of the Foldy-Twersky equation can be found [17,19] in terms of the complex forward-scattering amplitude of the aggregate. As regards the propagation of the coherent field, the fiber samples behaves as a homogeneous medium with a complex refractive index. Note that, even if the spherical particles that constitute the aggregate have a real refractive index and are therefore non-absorptive, the refractive index defined through this theoretical approach is complex, so that the dispersion as a whole behaves like an absorbing medium [17].
Thus, we can calculate the albedo at =510 nm for systems with increasing thickness. Results from Tmatrix calculations for each system are shown as open circles in Figure S10. Thicker samples yield an albedo increasing towards the saturation value of A=1, for which the light is totally back-scattered. Following the approach given by Mishchenko et al. [7], a simple model can be used to follow the albedo dependence on sample thickness, that considers an exponential extinction according to the Lambert-Beer relation. In fact, we can consider a sigmoid dependence A(z)=A0exp(z/l*)/[1-A0+A0exp(z/l*)], where the intercept at zero thickness, A0=0.668, is fixed so that the value for the single layer (with thickness = 700 nm) matches the T-matrix calculations (A1=0.71) and l*= ls/(1-g)=3.4 μm (shown as a dashed line in Fig.  S10) is the transport length obtained from the similarity relation [19] presented in the main text. The model is shown as a blue line in Fig. S10, and it is in good agreement with the T-matrix calculations (open circles). One can conclude that light is highly extinct within few microns of propagation in the sample and for a thickness of 25 m the albedo is very close to unity.

Experimental albedo values.
The experimental albedo is obtained as the ratio of the normalized backscattered intensity (integrated over the backscattering angle) Iexp, and extinction, that is expressed in terms of transmission, Eexp=1-Tr=0.996. The integrated backscattered intensity was measured through integrated sphere measurements at =510 nm (Iexp=0.773). This yields a measured albedo, Aexp=0.8, for samples with thickness  25 m, that also approaches unity though being smaller than the corresponding calculated value (blue line in Fig. S10). This discrepancy can be justified by the experimental non-woven size which is much larger than the x-y size of modelled structure (10×10 µm 2 ), and by the fact that during measurements a part (3%) of light is lost in the backward direction, since the used integrating sphere does not collect photons in a solid angle of 1.1×10 -3 sr around the direction of the incident beam (see Methods).