## Sample preparation

The marine diatoms *Skeletonema pseudocostatum* are isolated with the capillary pipette method33 from the Sarno River mouth (40.72875N, 14.466432E, Naples, Italy), and identified by 18s and 28s rDNA genes. Microalgae (initial concentration: 900000 cells/mL) are maintained in a f/2 medium34 prepared with sterile filtered (0.22 µm) natural seawater collected from oligotrophic areas of the Gulf of Naples. Diatoms are exposed to 0 (‘control’ sample), 10 or 100 µM of copper. In the latter cases, copper concentration is obtained by dissolving appropriate amounts of copper sulfate pentahydrate to the culture medium.

The concentration of living diatoms in all samples is estimated by cell counting with a Bürker chamber32. Each experiment is performed in triplicate. 10 ml of each sample are taken immediately after adding metals in the culture medium for measuring the THz absorption spectrum as explained in the next subsections. A 10 ml aliquot of the culture medium is employed as baseline for THz measurements.

## Tds-thz Experiment

The light pulses used for THz generation are provided by a regeneratively amplified femtosecond laser (Coherent Legend), which delivers optical pulses of 35 fs duration. The central wavelength of the laser is ∼800 nm and the bandwidth (full width at half-maximum) is ∼80 nm. The main part of the laser beam passes through the beam-splitter shown in Fig. 4a, where a diagram of our THz spectrometer is shown, and is focused by a lens. After the lens the beam passes through a β-barium borate nonlinear crystal that converts part of the fundamental beam in its second harmonic. In the focal point these intense beams generate a plasma formed by ionized air. Here, broadband and ultrashort THz pulses are generated via four-wave-mixing of the fundamental and second-harmonic photons35. The THz bandwidth is limited only by the laser bandwidth and therefore the highest generated frequency is 40 THz. The much less intense beam reflected by the beam-splitter is sent to an adjustable delay line and used for probing in time the generated THz pulse by means of the electro-optic effect in a suitable nonlinear crystal21. The latter technique allows a coherent detection of the THz electric field (amplitude and phase). Despite our very large THz bandwidth, the actual bandwidth available for the experiment is limited by the choice of the detection crystal. In the present work, we use a nonstandard crystal to reach efficient detection up to 10 THz: LAPC36. The spectrometer is closed in a chamber filled with high-purity (99.999%) nitrogen to strongly reduce the signal absorption by humid air. In Fig. 4b an example of power spectrum of a THz pulse propagating in nitrogen, transmitted through the most absorbing sample, and detected by electro-optic sampling in LAPC, is reported in logarithm scale. In Fig. 4c the corresponding time waveform of the optical field is shown. Despite the bandwidth of our detection crystal, the sample absorption drastically reduces the transmitted signal, so that the signal-to-noise ratio for the experiments is good only up to about 5 THz in the most absorbing samples. Yet, this bandwidth is not usual for THz-TDS applied to aqueous solutions. Liquid samples are placed in a self-assembled variable-path cell with silicon access windows. Silicon is used because it has a good and flat transmission up to 20 THz. The cell path may be continuously varied from 1 mm down to few microns with an average step of about 14 \(\mu\)m. For technical reasons, the ruler ticks are not equally spaced. However, each step is finely calibrated by measuring the temporal delay from the main pulse of the THz pulse replica transmitted after reflections from the two air-silicon interfaces internal to the empty cell. Then, we measure the THz pulses transmitted through the liquid samples for different thicknesses. For all measurements we use the following thicknesses: 22, 29, 47, 64, 79, and 94 \(\mu\)m.

**Absorption coefficient and refractive index extraction**

The THz fields \(E\left(t\right)\) measured for different thicknesses are transformed in the frequency domain by Fast Fourier Transform to obtain the complex field \(\tilde{E}\left(\nu \right)=E\left(\nu \right){e}^{i\varphi \left(\nu \right)}\), where \(\nu\) represents the frequency, \(E\left(\nu \right)\) and \(\varphi \left(\nu \right)\) represent the module and the phase of the complex field, respectively. As we will see in the following, these two quantities are directly linked to the complex refractive index, \(\tilde{n}\left(\nu \right)=n\left(\nu \right)+ik\left(\nu \right)\), of the sample. Here, \(n\) is the refractive index and \(k\) is the extinction coefficient. The latter is related to the absorption coefficient by \(\alpha =4\pi \nu k/c\), with \(c\) the speed of light in vacuum. In principle, the ratio of the Fourier transforms of only two THz temporal traces with and without the sample, with the second trace used as reference, provides the frequency-dependent complex refractive index \(\tilde{n}\left(\nu \right)\). However, the extraction of the complex dielectric function of liquid samples is complicated by the presence of the container, whose windows affect the THz signal as well. In this case, the complex refractive index is more reliably extracted by isolating the sample response from the window response through measurements at varying sample thickness \(d\). With this approach, the reference is the sample itself and the contribution of the cell is directly taken into account. In detail, the complex refractive index may be obtained by using the following Beer-Lambert generalized relationships:26

$$n\left(\nu \right)=c \frac{\varphi \left(\nu ,d+\varDelta d\right)-\varphi (\nu ,d)}{2\pi \nu \varDelta d}=\frac{\text{c}{\Delta }\varphi \left(\nu \right)}{2\pi \nu \varDelta d}$$

1

,

$$k\left(\nu \right)=c \frac{\text{ln}E\left(\nu ,d\right)-\text{ln}E(\nu ,d+\varDelta d)}{2\pi \nu \varDelta d}=\frac{c{\Delta }\text{ln}E\left(\nu \right)}{2\pi \nu \varDelta d}$$

2

.

For all our measurements *d* = 22 \(\mu\)m, i. e., the thinnest sample, while \(\varDelta d\) is the change in thickness of all the other samples. As it is apparent from Eqs. 1 and 2, the phase change \({\Delta }\varphi \left(\nu \right)\) and the change of the natural logarithm of the electric field module \({\Delta }\text{ln}E\left(\nu \right)\) depend linearly on \(\varDelta d\) for a given frequency. Therefore, at each frequency, a linear model can be used for fitting the experimental results to extract \(n\left(\nu \right)\) and \(k\left(\nu \right)\), which become the slopes of these linear models. To fit the latter we use a standard ‘chi-square’ minimization with the axis intercept of the linear model constrained to zero37.

If we assume that the statistical errors are the same for all the data for a given frequency, then the best estimates for the slope, \(m\), and its error, \({\sigma }_{m}\), are given by the following relationships:

$$m=\frac{\sum _{i}^{N}{{\varDelta d}_{i}{\Delta }\text{ln}E\left(\nu \right)}_{i}}{\sum _{i}^{N}{\left({\varDelta d}_{i}\right)}^{2}}$$

3

$${\sigma }_{m}=\sqrt{\frac{\sum _{i}^{N}{\left({{m\varDelta d}_{i}-{\Delta }\text{ln}E\left(\nu \right)}_{i}\right)}^{2}}{(N-1)\sum _{i}^{N}\varDelta {d}_{i}^{2}}}$$

4

,

where \(N\) is the number of data. The statistical errors given by Eq. 4 are represented in Figs. 2 and 3 as shaded areas (66% confidence level). From Eq. 4 it is evident that the statistical error depends on the frequency \(\nu\). This dependence may lead to the strong modulations of the shaded areas graphed in Figs. 2 and 3.

In conclusion, we note that equations 1 and 2 are valid only if multiple reflections in the sample thin layer are neglected (those occurring in the silicon windows are not considered since they appear out of the measurement temporal window). In case of thin samples this is possible only for very absorbing samples, as in our case, since the pulse is strongly attenuated during the back-and-forth trip so that the echoes of the main pulse become negligible.