# Waveguide QED at the onset of spin-spin correlations. Supplemental Material

Sebastián Roca-Jerat  Instituto de Nanociencia y Materiales de Aragón (INMA),
CSIC-Universidad de Zaragoza, Zaragoza 50009, Spain  Marcos Rubín-Osanz
Instituto de Nanociencia y Materiales de Aragón (INMA), CSIC-Universidad de
Zaragoza, Zaragoza 50009, Spain  Mark D. Jenkins  Instituto de Nanociencia y
Materiales de Aragón (INMA), CSIC-Universidad de Zaragoza, Zaragoza 50009,
Spain  Agustín Camón  Instituto de Nanociencia y Materiales de Aragón (INMA),
CSIC-Universidad de Zaragoza, Zaragoza 50009, Spain  Pablo J. Alonso
Instituto de Nanociencia y Materiales de Aragón (INMA), CSIC-Universidad de
Zaragoza, Zaragoza 50009, Spain  David Zueco  Instituto de Nanociencia y
Materiales de Aragón (INMA), CSIC-Universidad de Zaragoza, Zaragoza 50009,
Spain  Fernando Luis  Instituto de Nanociencia y Materiales de Aragón (INMA),
CSIC-Universidad de Zaragoza, Zaragoza 50009, Spain<EMAIL_ADDRESS>
##  I Sample preparation and structural characterization

DPPH ($2,2-$diphenyl$-1-$picrylhydrazyl) is a free radical molecule hosting a
spin $S=1/2$ with a nearly isotropic $g_{S}=2.004$ factor, very close to that
of a free electron.Williams (1966); Žilić  _et al._ (2010) This molecule is a
widely used standard in Electron Paramagnetic Resonance (EPR).Yordanov (1996)
DPPH samples used throughout this work were purchased, in powder form, from
Sigma Aldrich (reference D9132).

Even though solvent-free DPPH has been known for a long time,Williams (1966) a
full structural determination was not performed until quite recently.Žilić
_et al._ (2010) The results of powder X-ray diffraction experiments performed
on our sample are shown in Fig. S1. They suggest that it corresponds to the
DPPH-III structure,Žilić  _et al._ (2010) shown in Fig. 1a in the main text,
which contains two inequivalent DPPH sites, referred to as A and B
sublattices.

Figure S1: Determination of the crystal structure of the DPPH samples with
X-ray diffraction experiments. Powder diffraction experiment performed on a
DPPH sample (blue solid line) and the calculated pattern for the DPPH-III
crystal structure (red dashed line).Žilić  _et al._ (2010)

##  II Magnetic characterization

Figure S2: Micro-Hall magnetometry. (a) Optical microscopy image showing a
bare micro-Hall chip with two Hall crosses (left) and a chip hosting a DPPH
sample on top of one of them (right). (b) Schematic image of the magnetometer,
illustrating the differential detection method used to measure the Hall
voltage associated with the sample’s magnetization.

###  II.1 Magnetic properties measurements

The dc magnetic susceptibility and magnetic isotherms were measured above $2$
K in a Magnetic Properties Measurement System (MPMS) by Quantum Design, a
commercial magnetometer based on a dc-SQUID, operated by the Servicio de Apoyo
a la Investigación (SAI) of the University of Zaragoza. In addition,
measurements were also performed, for $0.3$ K $\leq T\leq 10$ K, with a
homemade micro-Hall magnetometer Kent  _et al._ (1994) mounted in the insert
of a Physical Properties Measurement System (PPMS), also by Quantum Design and
operated by the SAI. The magnetometer (Fig. S2) consists of two semiconducting
layers of GaAs and Al1-αGaαAs, shaped in the form of a double cross, with a
two-dimensional electron gas confined at the interface between the two
materials. A current $I_{\rm ac}$ is introduced through the main arm and a
magnetic field $\vec{B}$ is applied parallel to it. The DPPH sample was placed
on the edge of one of the two Hall crosses of the device and fixed/thermalized
with Apiezon N grease. Its magnetization $\vec{M}$ generates an additional
magnetic field $\vec{B}_{M}$, perpendicular to the plane of the device, which
induces a Hall voltage $V_{\rm Hall}$ proportional to $|\vec{M}|$. In order to
reduce noise and detect this signal, the current was modulated at a frequency
of $107$ Hz, with an $I_{\rm ac}=10$ $\mu$A amplitude, and the reference
voltage $V_{\rm ref}$ in the bare cross was subtracted from $V_{\rm Hall}$ by
means of a differential lock-in amplifier.

The dc-magnetic susceptibility $\chi$ measured between $0.35$ and $100$ K with
a magnetic field $B=0.1$ T can be found in Fig. 1 of the main text, as the
product $\chi T$. Figure S3 shows two complementary representations of these
data: $\chi$ and $1/\chi$ (Curie-Weiss plot) vs temperature. From room
temperature down to about $35$ K, $\chi$ follows a Curie-Weiss law
$\chi=C/(T-\theta)$, with $C$ the Curie constant and $\theta$ the Weiss
temperature. We find that $C=0.320(1)$ emu K/mol Oe measured in this region is
smaller, by a factor $x\equiv C/C_{S=1/2}\cong 0.85$, than $C_{S=1/2}=N_{\rm
A}g_{S}^{2}/4k_{\rm B}=0.375$ emu K/mol Oe that would be expected for isolated
$S=1/2$ spins having the g-factor $g_{S}=2.004$ reported for DPPH.Žilić  _et
al._ (2010) This suggests that a fraction of DPPH molecules are in their
oxidized form, which is diamagnetic, i.e. has $S=0$. The finite negative
$\theta=-21$ K reflects the existence of relatively strong antiferromagnetic
(AF) interactions. These interactions lead to the drop in the $\chi T$
product, proportional to the effective magnetic moment squared, that is
observed experimentally. A quite similar behaviour has been observed in
diverse DPPH derivatives.Ohya-Nishiguchi (1979); Fujito (1981); Žilić  _et
al._ (2010) It is associated with the coupling between nearest neighbour
molecules, which tend to form AF dimers with an $S=0$ ground state. In the
DPPH-III structure, such spin dimers are formed by DPPH molecules belonging to
the A crystal sublattice, see Fig. 1a in the main text and Žilić  _et al._ ,
2010.

Below $10$ K, $\chi\,T$ shows a second drop. In this temperature region,
$\chi$ follows also a Curie-Weiss law, but with a smaller $C\equiv C_{\rm
B}=0.154(1)$ emuK/mol Oe and a much smaller $\theta=-0.6$ K. This behaviour
probably reflects the weaker interactions between spins in sublattice B of
DPPH-III. The response remains paramagnetic down to very low-$T$, with no
clear indication of a phase transition. It is consistent with the formation of
either AF dimers or low dimensional ($1D$) antiferromagnetic chains. However,
specific heat data, discussed in section II.2, point to the formation of
chains. We have calculated the dc susceptibility of finite chains of $N\leq 8$
$S=1/2$ spins coupled by a Heisenberg AF coupling. Details of the model and
the calculations are given in section V.1 below. The experimental $\chi T$ can
be described by the combination of two contributions: one arising from AF
dimers of A-type DPPH molecules and another one resulting from AF chains
formed within the B sub-lattice. The reason why the low-$T$ $\chi$ can be
accounted for with a model that considers small-$N$ chains is that the
contribution per spin to the total susceptibility tends to vanish as
$N\rightarrow\infty$. The uncompensated spins in odd-$N$ chains are
responsible for the incomplete cancellation of $\chi T$ as $T\rightarrow 0$.
The origin of finite size chains can be associated with the presence of about
$15$% of diamagnetic DPPH molecules, which are randomly distributed over the
crystal lattice.

Magnetic isotherms measured at $0.5$, $1$, $2$, $5$ and $10$ K for magnetic
fields up to $8$ T are also shown in Fig. S3. Below 5 K, the theory for non-
interacting paramagnetic spins (solid lines) overestimates the observed
magnetization (light dots), which is partially suppressed at low fields by AF
interactions.

Figure S3: Magnetic characterization. (a) dc-magnetic susceptibility $\chi$ of
a DPPH powdered sample measured between $0.35$ and $300$ K with a magnetic
field of $0.1$ T. Open dots were measured with a dc-SQUID magnetometer, solid
dots with a micro-Hall magnetometer. The predictions of a model of infinite AF
chains (red solid line) and a model of AF dimers (cyan solid line) are
compared to the experimental data. The black dashed line is the Curie law for
a fully paramagnetic DPPH sample. (b) Curie-Weiss plot of the reciprocal
susceptibility $\chi^{-1}$ vs temperature. The light red dashed line is a
least square linear fit of data measured above $50$ K. It gives a Weiss
temperature $\theta=-21$ K associated with the formation of AF dimers in
sublattice A. Inset: close up plot of the low temperature region data, which
follow a Curie-Weiss law (light red dashed line) with $\theta=-0.65$ K
associated to the formation of AF chains in sublattice B. (c) Magnetic
isotherms, measured at 0.5, 1, 2, 5 and 10 K for magnetic fields up to 8 T..
The lines show the dependence expected for non-interacting paramagnetic spins.

###  II.2 Heat capacity measurements

Heat capacity experiments were carried out also in a PPMS. The calorimeter
consists of a sapphire disc — the sample holder— on which the sample is
placed. This holder integrates a heater and thermometer, and it is connected
to a thermal reservoir by thin gold wires. Sample and holder are put in mutual
thermal contact with apiezon N grease. Measurements are performed with
controlled pressure conditions in the sample chamber of the PPMS, that is,
constant-pressure heat capacity is measured. The heat capacity of the sample
is measured using a relaxation technique.Stewart (1983) A heat power is
applied for a short time, increasing the temperature of both the holder and
the sample. Part of the heat is transferred through the wires to the thermal
reservoir. The time constant of the exponential increase and decrease of
temperature during this process is proportional to the heat capacity of the
combined sample plus holder system. The different holders are regularly
calibrated without sample, so that the heat capacity of the holder can be
subtracted from the measurement.

Figure S4: Heat capacity of polycrystalline DPPH for different magnetic
fields. The experimental results (dots) are compared with a theory of small AF
chains (solid lines) and with the paramagnetic theory for non-interacting
$S=1/2$ spins (light dashed lines).

The zero-field specific heat of DPPH is shown in Fig. 1c of the main text.
Above $2$ K, it is dominated by a large contribution arising from lattice and
molecular vibrations. At lower temperatures, an additional contribution shows
up, which depends on magnetic field, as shown in Fig. S4. Therefore, it
reflects spin excitations. The broad shape of this anomaly confirms that no
phase transition to long-range magnetic order takes place in this temperature
region.

The zero-field data agree well with predictions for the specific heat of
infinite AF chains,Blöte (1975) and seem to rule out the formation of dimers
within the B sublattice. Lacking an analytical model for the $B\neq 0$
specific heat of infinitely long AF spin chains, we had to content ourselves
with computing the heat capacity of finite chains of up to eight $S=1/2$
spins. The model is discussed in section V.1. This simplified model accounts
reasonably well for data measured at low ($B\leq 0.5$ T and high ($B>2$ T)
magnetic fields. In the latter case, it tends towards the predictions for non-
interacting spins, as expected. Deviations are observed at intermediate fields
(1–2 T).

##  III Electron paramagnetic resonance experiments

Electron Paramagnetic Resonance (EPR) experiments were performed in a
commercial Elexsys E-580 by Brucker operating in the X-band ($\sim$9.8 GHz).
The spectrometer consists of a resonant cavity, in which the paramagnetic
sample is placed, located at the center of an electromagnet. The external DC
magnetic field generated by the magnet tunes the spin transition to the cavity
frequency, which results in a resonant absorption of microwave power by the
sample.

In continuous wave EPR (cw-EPR) experiments, the cavity is continuously
irradiated while the DC magnetic field is swept slowly. The absorption signal
is obtained with a field-modulation detection, in which a set of small coils
apply a radio-frequency (RF) magnetic field modulation that results in a
derivative-like line shape for each absorption signal. The single narrow line
shape measured with this technique, shown in Fig. S5a, agrees with a spin
$S=1/2$ and $g_{S}=2.004$. The line intensity $I$ shows a paramagnetic
behaviour, with $IT$, also shown in Fig. S5, behaving very much as $\chi T$.
The line shape suggests that the line width is dominated by the homogeneous
broadening down to $4$ K. It is known that, in the paramagnetic phase of DPPH,
exchange interactions tend to suppress the inhomogenous broadening associated
with dipole-dipole interactions Vleck (1948) via the exchange narrowing
mechanism.Anderson and Weiss (1953) This property makes DPPH a good reference
material for calibrating commercial EPR systems.Yordanov (1996)

Figure S5: cw-EPR experiments. (a) cw-EPR spectra of a powder DPPH sample
measured at five different temperatures. (b) Evolution of the product of the
intensity of the cw-EPR spectrum times temperature, measured first for
decreasing temperatures (blue closed dots) and then for increasing
temperatures (red open dots).

Time-domain, or pulsed, EPR experiments measure the signal emitted by the spin
system after sending a series of microwave pulses to the cavity. This gives
information on the relaxation and coherence of the spin system, as the
detected signal depends on the dynamics of the spin state prepared by the
pulses. The simplest pulsed EPR experiment is the measurement of the Free
Induction Decay (FID), which is the signal that the spin system induces in the
cavity after a microwave pulse generates a coherence between the two spin
states of a given transition. If the spin system exhibits inhomogeneous
broadening, a revival of the FID signal — the ‘echo’ — can be recovered with
refocusing pulses. When homogeneous broadening is dominant, there is no spin
echo signal, and the phase memory time $T_{\rm{m}}$ of the spin system is
obtained as the characteristic time constant of the FID. Time-domain EPR
measurements of the FID of DPPH powder samples were carried out at $7.7$ K and
at room temperature. The results are shown in Fig. S6. The fit of the two
signals is compatible with a homogeneous broadening $\gamma_{\phi}/2\pi=1/2\pi
T_{\rm{m}}=4.8$ MHz and $4.3$ MHz at 7.7 K and room temperature, respectively.
These values correspond to $T_{\rm{m}}=33$ ns and $37$ ns, respectively.
Figure S6c shows that these values are also compatible with those obtained
from microwave transmission experiments at lower temperatures.

Figure S6: Pulsed-EPR experiments. (a) Free Induction Decay (FID) of a DPPH
powered sample measured at 7.7 K. (b) FID of the same sample measured at room
temperature. (c) Temperature dependence of the DPPH absorption spectrum width.
The values extracted from FID data (red open dots) are compared to those
derived from on-chip microwave transmission experiments (blue closed dots).

##  IV On-chip microwave transmission experiments

###  IV.1 Experimental setup

Microwave transmission experiments were performed with on-chip superconducting
coplanar waveguides fabricated by optical lithography on a $100$ nm Nb thin
film deposited onto crystalline sapphire wafers. A $1$ mm wide DPPH pellet was
fixed on the $400$ $\mu$m wide central line by means of apiezon N grease, see
Fig. S7a. Figures S7b-c show how the superconducting chip is integrated in the
3He-4He dilution refrigerator that gives access to temperatures as low as $7$
mK. A cold finger (Fig. S7b) places the chip at the center of a $1$ T axial
superconducting magnet (Fig. S7c). The magnetic field is applied along the
$\hat{y}$ axis of the laboratory reference frame, which is parallel to the
wave guide. The device is thermalized to the mixing chamber of the dilution
refrigerator (Fig. S7d). The transmission $S_{21}$ through the device was
measured in the frequency range $10$ MHz $\leq\omega/2\pi\leq 14$ GHz with a
vector network analyzer. The input signal was attenuated by 50 dB before
reaching the chip, while the transmitted signal was amplified by 30 dB at 4 K
with a low-noise cryogenic amplifier.

Figure S7: Experimental setup for on-chip transmission experiments. (a)
Microscopy image of a DPPH powder sample on top of the waveguide. (b) Cold
finger hosting the superconducting chip. (c) Superconducting magnet in which
the cold finger is inserted, generating a magnetic field parallel to the
waveguide. (d) Scheme of the 3He-4He dilution refrigerator with a base
temperature $\simeq 7$ mK.

###  IV.2 Normalization of microwave transmission data

Figure S8: Normalization of transmission data. (a) Raw microwave transmission
data $S_{21}$ measured on a waveguide coupled to a powdered DPPH sample at
five magnetic field values. They are chosen in such a way that spin resonances
at $\Omega=\mu_{\rm B}g_{S}B/\hbar$ do not overlap. Some of these resonances
are masked by spurious modes in $S_{21,\,\rm{modes}}$. (b) Normalization of
the transmission data with Eq. (3) at the same magnetic fields. Only the part
of $S_{21,\,\rm{WG+DPPH}}(B)\simeq 1+\Delta S_{21}(B)$ that is close to
$\Omega$ (solid lines) is used to extract the spin-photon coupling $G$ and the
dissipation rate $\Gamma$. The rest (dashed lines), which includes an upwards
peak associated with the subtraction of $S_{21}(B+\Delta B)$ (in this example,
$\Delta B<0$), is discarded. This peak is not relevant as long as it does not
overlap with the absorption peak at $\Omega$, which sets the condition
$\mu_{\rm B}g_{S}|\Delta B|/h\gg\Gamma/2\pi$.

In the actual transmission data the microwave absorption by the DPPH spins is
often masked by spurious modes, see Fig. S8a, that arise from unavoidable
imperfections in the microwave circuit. The measured transmission can then be
understood as the convolution:

$S_{21}\simeq S_{21,\,\rm{WG+DPPH}}\,S_{21,\,\rm{modes}}\;,$ (1)

of the transmission features $S_{21,\,\rm{WG+DPPH}}$ arising from the presence
of the DPPH sample on top of the waveguide with the contribution
$S_{21,\,\rm{modes}}$ arising from those modes.

In order to obtain $S_{21,\,\rm{WG+DPPH}}$, the raw $S_{21}$ needs to be
normalized. A simple normalization procedure exploits the magnetic field
dependence of $S_{21,\,\rm{WG+DPPH}}$:

$S_{21,\,\rm{WG+DPPH}}(B)\simeq 1+\Delta S_{21}(B)\;,$ (2)

with:

$\Delta S_{21}(B)=\dfrac{S_{21}(B)-S_{21}(B+\Delta B)}{S_{21}(B_{\rm
ref})}\;.$ (3)

Here, the raw data $S_{21}(B)$ measured at a field $B$ are normalized by
subtracting from it the raw data $S_{21}(B+\Delta B)$ measured at a close
field $B+\Delta B$ and dividing the result by a reference transmission
measured at a fixed field $B_{\rm res}$ (usually $B_{\rm res}=0$ is chosen, or
any magnetic field at which the microwave transmission shows no features
associated to any spin resonances in the relevant frequency region). As shown
in Fig. S8b, this method effectively removes $S_{21,\,\rm{modes}}$ from the
transmission data, provided that $S_{21,\,\rm{modes}}$ is independent of
magnetic field or that it varies sufficiently slowly with it. The
normalization step $\Delta B$ is chosen to be as small as possible while
avoiding the overlap of the spin resonances in $S_{21}(B)$ and
$S_{21}(B+\Delta B)$. This requires that $|\Delta B|\gg\hbar\Gamma/\mu_{\rm
B}g_{S}$.

###  IV.3 Transmission properties in the paramagnetic phase of DPPH

Figure S9 shows the field and temperature dependence of $G$ and $\Gamma$ in
the paramagnetic phase of the B-type DPPH sublattice ($T>$ 0.7 K). The
coupling $G$ is proportional the input frequency in a one-dimensional wave
guide. Close to the spin resonance, this input frequency is approximately the
spin resonance frequency $\Omega=\mu_{\rm B}g_{S}B/\hbar$, and $G$ can be
written as (see Eq. (3) in the main text):

$G=2\pi\alpha\Omega N\tanh{\left(\dfrac{\hbar\Omega}{2k_{\rm B}T}\right)}\;,$
(4)

where $N$ is the total number of spins coupled to the microwave magnetic field
generated by the waveguide, and $\alpha$ is a temperature-independent and
field-independent dimensionless scaling factor. The only free parameter in the
fit of $G$ for all temperatures and magnetic fields within the paramagnetic
phase is $2\pi\alpha N=$ 0.00441(6). The coupling enhancement by
$N_{\rm{eff}}=N\tanh{\left(\hbar\Omega/2k_{\rm B}T\right)}$ signifies the
increase in the radiative rate decay in superradiant collective states.

Figure S9: Magnetic field dependence of $G$ and $\Gamma$ at four selected
temperatures within the paramagnetic phase of DPPH. $G$ (red open dots) and
$\Gamma$ (blue closed dots) values were obtained by fitting the normalized
transmission data using Eq. (2) in the main text. Light solid lines show the
average $\Gamma$ (discarding the low magnetic field region) and the result of
a least-squares fit of $G$ with Eq. (4).

For $T>1$ K, the dissipation rate $\Gamma$ is found to depend weakly on both
magnetic field and temperature. The values derived from these experiments are
compatible with those obtained with a conventional EPR spectrometer, see Fig.
S6c. Besides, $\Gamma/2\pi\cong 8$ MHz measured at $T=4.2$ K is also
compatible with that obtained for DPPH samples coupled to coplanar resonators
at the same temperature.Gimeno  _et al._ (2020) Remarkably, $G$ becomes
comparable to $\Gamma$ for magnetic fields $B>$ 400 mT. This yields a signal
visibility ${\rm max}(|\Delta S_{21}|)=G/(G+\Gamma)$ close to unity.
Therefore, the spin ensemble and the waveguide enter then a regime that
resembles the strong coupling limit ($G\geq\Gamma$) that is commonly reached
only by coupling spins to a resonator. In order to better understand how large
this $G$ actually is, it is worth estimating what it would lead to in the case
of a resonator. Input-output theoryHümmer  _et al._ (2013) provides a relation
between $G$ and the coupling $G_{\rm res}$ to a cavity resonating at a given
frequency $\omega$

$G_{\rm res}=\sqrt{\frac{G\omega}{\pi}}$ (5)

Using this relation, the maximum coupling $G/2\pi\simeq 12$ MHz to the line
that we measure at $B=500$ mT and $\omega/2\pi=14$ GHz would correspond to
$G_{\rm res}\sim 0.58$ GHz. This estimate compares favourably to those
reported previously for DPPH samples at coplanar superconducting
resonators.Ghirri  _et al._ (2016); Mergenthaler  _et al._ (2017)

###  IV.4 Effect of input microwave power

The resonance of the DPPH ensemble was measured at $T$ = 0.3 K and $B$ = 0.22
T with different input microwave powers $P_{\rm in}$ between -25 and -10 dBm.
No saturation effect was observed in the resonance due to the change in
$P_{\rm in}$, as shown in Fig. S10.

Figure S10: Spin resonance versus microwave power. Normalized transmission of
a waveguide coupled to a powdered DPPH sample measured at $T=0.3$ K and
$B=0.22$ T.

##  V Theory of antiferromagnetic spin chains

###  V.1 Exact diagonalization calculation of the magnetic and thermal
responses of small AF spin chains

The 1D antiferromagnetic chains that form in the B-sublattice can described by
($\hbar=1$):

$\mathcal{H}_{m}=\sum_{i,j}\vec{S}_{i}\cdot\bar{J}\cdot\vec{S}_{j}-g\mu_{B}B\sum_{i}S_{i}^{y}\;.$
(6)

Here, $\vec{S_{j}}=(\sigma_{j}^{x},\sigma_{j}^{y},\sigma_{j}^{z})$ are the
Pauli matrices. The anisotropy of the spin-spin interactions is encoded in the
coupling tensor $\bar{J}$:

$\bar{J}={\rm diag}(J,J(1+\epsilon\sin\psi),J(1+\epsilon\cos\psi))\;.$ (7)

Here, $\psi$ is the angle that the anisotropy axis $\vec{\epsilon}$ forms with
$z$. Moreover, w.l.o.g we assume this axis lies in the $yz$-plane, see figure
S11. The last term in $\mathcal{H}_{m}$ is the Zeeman coupling to the
externally applied DC magnetic field $B_{0}$. We define the $y$-axis as the
orientation of this field, matching the laboratory frame defined in Fig. 1a in
the main text.

Figure S11: Orientations of the magnetic field and of the magnetic anisotropy
axis, and angle convention used in the exact diagonalization of the
Hamiltonain of small 1D AF chains.

Using exact diagonalization we can obtain all the eigenstates of the system up
to sizes not much larger than 10 spins. With them, we can calculate
observables such as the specific heat or the magnetic susceptibility. Starting
with the specific heat, $c$, we calculate it from the eigenvalues of the
Hamiltonian given in Eq. (6) at a certain temperature, $T$, with the following
expression

$c=\beta\left[\frac{Z^{\prime\prime}}{Z}-\left(\frac{Z^{\prime}}{Z}\right)^{2}\right]\
,$ (8)

where $\beta=1/T$, $()^{\prime}=\partial_{\beta}$ and
$Z=Tr\left(e^{-\beta\mathcal{H}}\right)$.Garcia-Palacios  _et al._ (2009) The
terms explicitly read $Z=\sum_{m}e^{-\beta\epsilon_{m}}$,
$Z^{\prime}=-\sum_{m}e^{-\beta\epsilon_{m}}\epsilon_{m}$ and
$Z^{\prime\prime}=\sum_{m}e^{-\beta\epsilon_{m}}\epsilon^{2}_{m}$, with
$\epsilon_{m}$ being the eigenstates of Eq. (6). We compute these energy
fluctuations for the different angles $\psi$ and sizes up to $N=8$ spins and
then we take the average. The results of the simulations are presented in
Figs. S4 and 1c (main text). We confine our analysis to small-sized lattices,
as the complete spectrum is necessary to compute the specific heat at finite
temperatures, as discussed in Fig. S12. This limits the quantitative approach
to actual experimental values.

Figure S12: Specific heat obtained through exact diagonalization. The results
shown are for a 7-spin chain, as a function of temperature fixing the external
magnetic field at $0.5$ T. As can be seen, we only recover the analytical
curve (black curve in both plots) when the totality of the system’s
eigenstates is considered.

On the other hand, for the magnetic susceptibility, we can treat
perturbatively the partition function $Z$. In general, for a system described
by some hamiltonian $\mathcal{H}_{0}$ that is under a perturbation $Q$ with
intensity $f$, i.e. $\mathcal{H}=\mathcal{H}_{0}-fQ$, the susceptibility,
$\chi$, is defined as

$\chi=\beta\left[\frac{Z_{2}}{Z_{0}}-\left(\frac{Z_{1}}{Z_{0}}\right)^{2}\right],$
(9)

where $Z\simeq Z_{0}+\xi Z_{1}+\xi^{2}Z_{2}$ and $\xi=\beta f$.Garcia-Palacios
_et al._ (2009) In our case, $Z_{0}=\sum_{m}e^{-\beta\epsilon_{m}}$,
$Z_{1}=\sum_{m}e^{-\beta\epsilon_{m}}Q_{mm}$, and
$Z_{2}=\sum_{mn}e^{-\beta\epsilon_{mn}}|Q_{mn}|^{2}K_{mn}$, with
$Q=\tfrac{1}{|\vec{B}|}\vec{B}\vec{S}$,
$K_{mn}=K\left[\beta(\epsilon_{m}-\epsilon_{n})\right]$ and
$K(X)=(e^{X}-1)/X$.

In this case, we observe that small chains accurately reproduce the
experimental curve (see Fig. S3). This phenomenon can be attributed to the
fact that, with increasing chain length, the per-spin contribution to the
total susceptibility diminishes. As demonstrated in Fig. S13, the
susceptibility scales as $1/N$, in relation to the susceptibility of a single
spin, which predominantly influences the behavior in the limit $T\rightarrow
0$. Furthermore, only chains with an odd number of spins contribute, due to
their $S=1/2$ ground state. At sufficiently high temperatures, however, chains
of all lengths contribute equally, enabling the recovery of the experimental
results using solely small chains, as shown in Fig. 1b in the main text.

Figure S13: Dependence of the magnetic susceptibility on the spin chain size.
We show the $\chi T$ product calculated, at $B=0$, as a function of the number
of spins in the chain in two regimes: low temperature (blue squares) and high
temperature (red triangles). In the low-temperature regime only odd number
spin chains contribute. This response (normalized by the spin number) decays
with increasing chain length.

Comparing the specific heat and the magnetic susceptibility obtained from the
exact diagonalization of our model of small chains with the experimental
results we found that the anisotropy in the spin-spin interactions is very
small and not relevant in these experiments. This is not the case in microwave
transmission experiments, where the small anisotropy plays a crucial role and
significantly changes the spin resonance features, as it is shown in the main
text and in section VI below.

###  V.2 Mean field theory of antiferromagnetic spin chains

For computing the transmission at very low temperatures, it is necessary to
consider large chains that are beyond the scope of exact diagonalization.
Therefore, we employ a mean field theory, which not only explains the
underlying physics but also accounts for the experimental data both
qualitatively and quantitatively, as detailed in the main text. Here, we aim
to describe our application of this mean field (MF) theory to DPPH.
Additionally, within the MF theory, we can compute the linear response,
enabling us to characterize the transmission experiments.

We use the Bogolyubov inequality, which serves as a method for approximating
the equilibrium partition function of an interacting system. This inequality
arises from convexity and requires the Hamiltonian to be divided into a part
that is solvable and the remaining part.

${\cal H}={\cal H}_{0}+{\cal H}_{1}\;\Longrightarrow\;\frac{Z}{Z_{0}}\geq{\rm
e}^{-\beta\langle H_{1}\rangle_{0}}\;,$ (10)

with $Z={\rm Tr}[{\rm e}^{-\beta{\cal H}}]$ the actual partition function,
$Z_{0}={\rm Tr}[{\rm e}^{-\beta{\cal H}_{0}}]$ and $\langle...\rangle_{0}={\rm
Tr}[...{\rm e}^{-\beta{\cal H}_{0}}]$. Taking the logarithm in Eq. (10), an
upper bound for the free energy is found:

$F\leq F_{0}+\langle{\cal H}_{1}\rangle_{0}\;.$ (11)

In antiferromagnetic systems, spins are divided in two sub-lattices, namely
$1$ and $2$. In each sub-lattice the spin averages are the same. Therefore, we
choose for ${\cal H}_{0}$:

${\cal H}_{0}=-\vec{\lambda}_{1}\vec{M}_{1}-\vec{\lambda}_{2}\vec{M}_{2}\;,$
(12)

where $\vec{\lambda}_{\alpha}$ are variational parameters that minimize the
left hand side of Eq. (11) and we have used the following compact notation for
the magnetization vectors:

$\displaystyle\vec{M}_{1}=$ $\displaystyle M_{1}\vec{u}_{1}\;,$ (13)
$\displaystyle\vec{M}_{2}=$ $\displaystyle M_{2}\vec{u}_{2}\;.$ (14)

It turns out that $Z_{0}$ only depends on the modulus
$|\vec{\lambda}_{\alpha}|$. Thus, we can set $\lambda_{\alpha}\geq 0$.
Considering nearest neighbor interactions, together with the angle convention
used in Fig. S11, we obtain for the average:

$\displaystyle\langle{\cal H}_{1}\rangle_{0}=\;$
$\displaystyle\lambda_{1}M_{1}+\lambda_{2}M_{2}$ (15)
$\displaystyle-g\mu_{\mathrm{B}}B(M_{1}\sin\theta_{1}+M_{2}\sin\theta_{2})$
$\displaystyle+M_{1}M_{2}\vec{u}_{1}\cdot\bar{J}\cdot\vec{u}_{2}\;,$

where $\bar{J}$ is given in Eq. (7). We assume that the magnetization vectors
lie also in the $zy$-plane, the plane spanned by the anisotropy and magnetic
field vectors as depicted in Fig. S11. Simulations performed considering a
generic three-dimensional orientation for the magnetisation vectors gave the
same results as in the 2D case. This is easy to see since the directions of
interest are always going to be the anisotropy axis (when there is no field)
or the field axis (when the field is very strong).

Inserting Eq. (15) in Eq. (11) and minimizing over the variables
$\\{\lambda_{1},\theta_{1},\lambda_{2},\theta_{2}\\}$, we end with the
consistency relations:

$\displaystyle\frac{\partial\langle{\cal H}\rangle_{0}}{\partial\theta_{1}}$
$\displaystyle=0\;,$ (16a) $\displaystyle\lambda_{1}$
$\displaystyle=g\mu_{\mathrm{B}}B\sin\theta_{1}-M_{2}\vec{u}_{1}\cdot\bar{J}\cdot\vec{u}_{2}\;,$
(16b) $\displaystyle M_{1}$ $\displaystyle=\tanh(\beta\lambda_{1})\;,$ (16c)

and similar equations for sublattice 2.

In general, these equations are solved numerically. However, analytical
progress can be made using symmetry arguments. We assume
$\lambda_{1}=\lambda_{2}=\lambda$ and $M_{1}=M_{2}=M$. With our convention for
angles (Fig. S11, we find that in either the paramagnetic phase or the spin-
flop transition (see below), $\theta_{1}+\theta_{2}=\pi$. This condition is
also applicable, in principle, to the antiferromagnetic scenario. By imposing
these conditions, we derive the following expression:

$\sin\theta_{2}^{\mathrm{eq}}=\frac{g\mu_{\mathrm{B}}B}{MJ[2+\epsilon(\sin\psi+\cos\psi)]}\;,$
(17)

with:

$\displaystyle M=\tanh\bigg{[}\beta\big{(}$ $\displaystyle
g\mu_{\mathrm{B}}B+JM\cos 2\theta_{2}^{\mathrm{eq}}$ (18) $\displaystyle-
JM\epsilon(\sin^{2}\theta_{2}^{\mathrm{eq}}\sin\psi-\cos^{2}\theta_{2}^{\mathrm{eq}}\cos\psi)\big{)}\bigg{]}\;.$

Equations (17) and (18) are used for determining the phase diagram, as
illustrated in Fig. S14. Here, we depict the numerical solutions for the
equilibrium angles, $\Delta\theta_{\mathrm{eq}}$, yielding the distinct
phases. A zero value of $\Delta\theta_{\mathrm{eq}}$ indicates that
temperature and/or magnetic field are sufficient to align the magnetization
vectors. A non-zero $\Delta\theta_{\mathrm{eq}}$ signals the emergence of
magnetic correlations, revealing two distinct phases: the spin-flop phase for
$0<\Delta\theta_{\mathrm{eq}}<\pi$, and the antiferromagnetic phase for
antiparallel spins. Panel (a) of Fig. S14 displays the case with perpendicular
anisotropy to the magnetic field, while panel (b) considers parallel
orientation. Specifically, when $\vec{B}\perp\vec{\epsilon}$ and $\psi<\pi/4$,
a magnetic field weaker than $J\sqrt{2\epsilon(\sin\psi-\cos\psi)}$ results in
antiparallel magnetization vectors orthogonal to the anisotropy axis, where
$\theta_{1}+\theta_{2}=2\pi$. This condition persists until the magnetic field
$B$ is strong enough to rotate both vectors towards the paramagnetic phase.
This process is elaborated in subplots (a.1) and (b.1) of Fig. S14. The
critical field at $T=0$ is derived from Eq. (17) as
$g\mu_{\mathrm{B}}B^{\mathrm{c}}=J(2+\epsilon(\sin\psi+\cos\psi))$, while the
Néel temperature is fixed at $T_{\mathrm{N}}=J/k_{\mathrm{B}}$.

Figure S14: Equilibrium magnetization angles. Phase diagram for
$\Delta\theta_{\rm eq}$ as a function of temperature and magnetic field. In
(a) the anisotropy is perpendicular to the magnetic field while in (b) they
are parallel. In subplots (a.1) and (b.1) we see the different behaviour at
$T=0$ that the anisotropy shows up. Since $J(1+\epsilon)<J$, the magnetization
vectors lie perpendicular to the anisotropy axis to minimize the energy. When
$\psi<\pi/4$, there is a magnetic field threshold for which, below that point,
magnetization remains in the perpendicular axis to $\vec{\epsilon}$. When the
field is strong enough, both vectors flip to start coming parallel between
them symmetrically until the paramagnetic phase is reached. This point is
found to be $g\mu_{\rm B}B=J\sqrt{2\epsilon(\sin\psi-\cos\psi)}$. The
parameters used were $J/k_{\rm B}=0.7\ \rm K$ and $J\epsilon/k_{\rm B}=-0.065\
\rm K$ for panels (a) and (b). In (a.1) and (b.1) an arbitrarily bigger
$\epsilon$ was used for better visualization.

###  V.3 Spin dynamics: computing the resonance frequency of the chains

Mean field treatment is consistent with using the Landau-Lifshitz-Gilbert
(LLG) equation for the dynamics. The LLG equation reads

$\frac{d\vec{M}_{\alpha}}{dt}=-g\mu_{B}\;\vec{M}_{\alpha}\times\vec{B}_{{\rm
mol},\alpha}-\gamma\;\vec{M}_{\alpha}\times(\vec{M}_{\alpha}\times\vec{B}_{{\rm
mol},\alpha})\,$ (19)

with $\alpha=1,2$. The molecular field is given by Vonsovskii (1966); Ross
(2013)

$\vec{B}_{{\rm mol},\alpha}=-\frac{1}{g\mu_{B}}\vec{\nabla}\,F\,,$ (20)

with $F$ the free energy computed by means of the MF, Cf. Eq. (11).

As the LLG preserves $|\vec{M}_{\alpha}|$, it is convenient to work in
spherical coordinates, which transform Eq. (19) into the following:

$\begin{split}\sin\theta_{\alpha}\;\dot{\theta}_{\alpha}&=-\frac{1}{M_{\alpha}}F_{\varphi_{\alpha}}-\gamma\sin\theta_{\alpha}\;F_{\theta_{\alpha}}\;;\\\
\sin\theta_{\alpha}\;\dot{\varphi}_{\alpha}&=\;\;\;\frac{1}{M_{\alpha}}F_{\theta_{\alpha}}-\gamma\frac{1}{\sin\theta_{\alpha}}F_{\varphi_{\alpha}}\;.\end{split}$
(21)

We have introduced the notation
$F_{\theta_{\alpha}}\equiv\partial_{\theta_{\alpha}}F$. At equilibrium, the
derivatives of the free energy are zero. Therefore, in the linear response
regime, we perform a Taylor expansion of these derivatives, resulting in the
set of linear equations:

$\displaystyle\begin{pmatrix}\sin\theta_{1}\;\dot{\theta}_{1}\\\
\sin\theta_{1}\;\dot{\varphi}_{1}\\\ \sin\theta_{2}\;\dot{\theta}_{2}\\\
\sin\theta_{2}\;\dot{\varphi}_{2}\\\ \end{pmatrix}=$
$\displaystyle\begin{pmatrix}-\frac{F_{\varphi_{1},\theta_{1}}}{M_{1}}&-\frac{F_{\varphi_{1},\varphi_{1}}}{M_{1}}&-\frac{F_{\varphi_{1},\theta_{2}}}{M_{1}}&-\frac{F_{\varphi_{1},\varphi_{2}}}{M_{1}}\\\
\frac{F_{\theta_{1},\theta_{1}}}{M_{1}}&\frac{F_{\theta_{1},\varphi_{1}}}{M_{1}}&\frac{F_{\theta_{1},\theta_{2}}}{M_{1}}&\frac{F_{\theta_{1},\varphi_{2}}}{M_{1}}\\\
-\frac{F_{\varphi_{2},\theta_{1}}}{M_{2}}&-\frac{F_{\varphi_{2},\varphi_{1}}}{M_{2}}&-\frac{F_{\varphi_{2},\theta_{2}}}{M_{2}}&-\frac{F_{\varphi_{2},\varphi_{2}}}{M_{2}}\\\
\frac{F_{\theta_{2},\theta_{1}}}{M_{2}}&\frac{F_{\theta_{2},\varphi_{1}}}{M_{2}}&\frac{F_{\theta_{2},\theta_{2}}}{M_{2}}&\frac{F_{\theta_{2},\varphi_{2}}}{M_{2}}\end{pmatrix}$
(22)
$\displaystyle-\gamma\begin{pmatrix}-\sin\theta_{1}\;F_{\theta_{1},\theta_{1}}&-\sin\theta_{1}\;F_{\theta_{1},\varphi_{1}}&-\sin\theta_{1}\;F_{\theta_{1},\theta_{2}}&-\sin\theta_{1}\;F_{\theta_{1},\varphi_{2}}\\\
\frac{F_{\varphi_{1},\theta_{1}}}{\sin\theta_{1}}&\frac{F_{\varphi_{1},\varphi_{1}}}{\sin\theta_{1}}&\frac{F_{\varphi_{1},\theta_{2}}}{\sin\theta_{1}}&\frac{F_{\varphi_{1},\varphi_{2}}}{\sin\theta_{1}}\\\
-\sin\theta_{2}\;F_{\theta_{1},\theta_{2}}&-\sin\theta_{2}\;F_{\theta_{2},\varphi_{1}}&-\sin\theta_{2}\;F_{\theta_{2}\theta_{2}}&-\sin\theta_{2}\;F_{\theta_{2},\varphi_{2}}\\\
\frac{F_{\varphi_{2},\theta_{1}}}{\sin\theta_{2}}&\frac{F_{\varphi_{1},\varphi_{2}}}{\sin\theta_{2}}&\frac{F_{\varphi_{2},\theta_{2}}}{\sin\theta_{2}}&\frac{F_{\varphi_{2},\varphi_{2}}}{\sin\theta_{2}}\end{pmatrix}\begin{pmatrix}\theta_{1}\\\
\varphi_{1}\\\ \theta_{2}\\\ \varphi_{2}\\\ \end{pmatrix}.$

Here, both $\theta_{\alpha}$ and $\varphi_{\alpha}$ are the equilibrium values
(obtained as explained in the previous section). Besides, the derivatives are
evaluated in the corresponding equilibrium values. Note that here we use a 3D
orientation for the magnetization vectors for convenience. The 2D choice that
we made in the previous section allowed us to derive some analytical
expression whereas here we will resort to numerical solutions.

The resonance frequency, $\Omega$ is derived by seeking oscillatory solutions
of the form:

$\displaystyle\theta_{\alpha}$
$\displaystyle=\delta\theta_{\alpha}\mathrm{e}^{i\Omega t},$ (23a)
$\displaystyle\varphi_{\alpha}$
$\displaystyle=\delta\varphi_{\alpha}\mathrm{e}^{i\Omega t}.$ (23b)

Substituting these expressions into Eq. (22), we obtain a linear equation
system for $\Omega$ that is amenable to numerical solution. Among the four
solutions, two correspond to the normal modes of each sublattice, indicating
that for each sublattice, two frequencies describe the motion of
$\theta_{\alpha}$ and $\varphi_{\alpha}$.

To calculate these resonant frequencies, we set $\gamma=0$, i.e., focusing
solely on Hamiltonian dynamics. While there is generally a dependence on the
dissipation strength, it is marginal. In Fig. (S15), we present several
numerical examples to illustrate the dependence of $\Omega$ on the angle
$\vec{\psi}$, which is crucial for understanding experimental outcomes.
Specifically, in panel S15a, we plot the angle dependence of the non-vanishing
resonant frequency as a function of the field for different anisotropy
orientations, $\psi$. The resonant frequencies are found to follow:

$\Omega=\sqrt{\bar{B}^{2}-2J^{2}\epsilon(\sin\psi-\cos\psi)},$ (24)

where $\bar{B}=g\mu_{\mathrm{B}}B$. For angles less than $\pi/4$, grey dashed
lines show the expected frequency curve absent a spin-flop transition.
Orientations with complementary angles ($\psi_{1}+\psi_{2}=\pi/2$) exhibit
identical resonant frequencies at zero field. The black dashed line denotes
the magnetic field value for which panel S15b is simulated, corresponding to
the experimental and simulation data presented in Fig. 3 in main text. Below
the Néel temperature ($T_{\mathrm{N}}$), magnetic correlations lead to
distinct resonant frequencies for each anisotropy orientation, resulting in a
distribution where larger angles, i.e., orientations closer to the magnetic
field direction, exhibit higher frequencies.

Additionally, the variation in resonance frequency due to magnetic
correlations necessitates estimating the line width (losses). Incorporating
the dissipation component, which is the second line in Eq. (22), and
considering the exponential decay in oscillatory solutions:

$\displaystyle\theta_{\alpha}$
$\displaystyle=\delta\theta_{\alpha}\mathrm{e}^{i\Omega t-\gamma t},$ (25a)
$\displaystyle\varphi_{\alpha}$
$\displaystyle=\delta\varphi_{\alpha}\mathrm{e}^{i\Omega t-\gamma t},$ (25b)

Moreover, numerical analyses suggest the approximation:

$\gamma\propto JM_{\alpha}^{2}\ .$ (26)

Thus, the line width increases below $T_{\mathrm{N}}$. Although this does not
determine the experimental broadening, which is due to the fact that
anisotropy introduces a distribution of resonance frequencies when spin waves
emerge, as explained in the main text.

Figure S15: Resonant frequency of the chains according to its orientation. In
a) we show the dependence of $\Omega$ at $T=0$ with the anistropy orientation,
$\psi$. We find that $\Omega$ follows Eq. (24) so angles between $0$ and
$\pi/4$ would follow the grey dashed line if they did not pass through the
particular spin flop transition described. The black dashed line indicates the
magnetic field at which the results in panel b) were obtained. There, the
temperature behaviour of the frequencies for each orientation $\psi$ is shown.
Over $T_{\mathrm{N}}$, the system is paramagnetic and anisotropy does not play
any role but below that temperature, magnetic correlations show up its strong
influence on the frequencies. The parameters used were $J/k_{\mathrm{B}}=0.7\
\mathrm{K}$ and $J\epsilon/k_{\mathrm{B}}=-0.065\ \mathrm{K}$.

##  VI Waveguide QED theory

The $B_{\rm rms}$ field is in the $xz$ plane, that is, perpendicular to the
applied field. Therefore, the photons from the guide will essentially produce
spin-flips in the DPPH. Thus, the Hamiltonian for our calculations will be
$\mathcal{H}={\mathcal{H}}_{\rm WG}+{\mathcal{H}}_{\rm I}+{\mathcal{H}}_{\rm
M}$, where

${\mathcal{H}}_{\rm WG}=\int{\rm
d}\omega\,\omega\,(r_{\omega}^{\dagger}r_{\omega}+l_{\omega}^{\dagger}l_{\omega})\
,$ (27)

${\mathcal{H}}_{\rm
I}=\sum_{j,n}^{N}\int_{\Omega_{-}}^{\Omega_{+}}d{\omega}\;\lambda_{\omega}({\vec{r}}_{j})\;\left(\sigma^{+}_{j,n}+{\rm
h.c.}\right)\,X(\omega)\ .$ (28)

In Eq. (27), $r_{\omega}^{\dagger}$ ( $l_{\omega}^{\dagger}$ ) are bosonic
operators,
$[r_{\omega},r^{\dagger}_{\omega^{\prime}}]=\delta(\omega-\omega^{\prime})$
and
$[l_{\omega},l^{\dagger}_{\omega^{\prime}}]=\delta(\omega-\omega^{\prime})$.
They create right (left) moving photons. Besides, in Eq. (28),
$X(\omega)=(r_{\omega}+l_{\omega})+{\rm h.c.}$, $\lambda_{\omega}$ are spin-
photon coupling constants and the $\sigma^{\pm}_{j,n}$ are Pauli spin-raising
and lowering operators for the $j$-th DPPH molecule from the $n$-th chain. The
angle $\theta$ stands for the orientation of the magnetic field $\vec{B}_{\rm
rms}$ generated by the waveguide at the position of the given molecule. The
actual value of the coupling is given by,

$\lambda_{\omega}(\vec{r}_{j})=\frac{g_{S}\mu_{\rm
B}}{\sqrt{2\pi}}|\vec{B}_{\rm rms}(\vec{r}_{j})|\cong\lambda_{\omega}$ (29)

with $\vec{B}_{\rm rms}(\vec{r}_{j})$ the zero-point root mean square magnetic
field generated by the current fluctuations in the waveguide. The last
approximation states that we neglect distance effects. Furthermore, the
dependence on the $x$ and $z$ coordinates can be averaged out to yield an
average coupling which, ultimately, must be fitted.

###  VI.1 Input-output fields and scattering

For computing the transmission, we use input-output theory Gardiner and
Collett (1985)

$r_{\rm out}(t)=r_{\rm in}(t)-i\int_{0}^{\infty}\frac{{\rm
d}\omega}{\sqrt{2\pi}}{\rm
e}^{-i\omega\,t}\lambda_{\omega}\int_{-\infty}^{\infty}{\rm d}\tau\;{\rm
e}^{+i\omega\tau}\sigma_{x}(\tau)$ (30)

With right _input_ photon fields

$\displaystyle r_{\rm in}(t)$ $\displaystyle:=\int_{0}^{\infty}\frac{{\rm
d}\omega}{\sqrt{2\pi}}\;r_{\omega}(t_{0}){\rm e}^{-i\omega(t-t_{0})}\,,$ (31)

and the _output_ ones,

$\displaystyle r_{\rm out}(t)$ $\displaystyle:=\int_{0}^{\infty}\frac{{\rm
d}\omega}{\sqrt{2\pi}}r_{\omega}(t_{f}){\rm e}^{-i\omega(t-t_{f})}$ (32)

and identical relations for the left fields.

####  VI.1.1 Input fields and temperature effects

We must consider the effect of temperature. The line-molecules are thermalized
at temperature $T$. Then, we _inject_ a current, that in the case of our
experiment is equivalent to a coherent state. Therefore the input field can be
written as,

$\varrho_{\rm in}={\mathcal{A}}\,\varrho_{\beta}\,{\mathcal{A}}^{\dagger}\;,$
(33)

with $\mathcal{A}$ the generator of coherent states,

$\displaystyle{\mathcal{A}}=\sum_{j}\phi_{j}{\rm
e}^{\alpha_{j}a_{j}^{\dagger}-\alpha_{j}^{*}a_{j}}\;.$ (34)

####  VI.1.2 Transmission and reflection formulas

In the limit $\alpha\to 0$ (cf. Eq. (34)), the transmission is given by:

$t^{(1)}(\omega)=\lim_{\alpha\to 0}\frac{\langle\alpha_{\omega}|r_{\rm
out}(t)|\alpha_{\omega}\rangle}{\langle\alpha_{\omega}|r_{\rm
in}(t)|\alpha_{\omega}\rangle},$ (35)

where $|\alpha_{\omega}\rangle=\exp{\left(\alpha r_{\rm
in}^{\dagger}(\omega)-\text{H.c.}\right)}|\Omega\rangle$, with $r_{\rm
in}^{\dagger}(\omega)$ being the Fourier transform of $r_{\rm
in}^{\dagger}(t)$. It is crucial to note that, in this limit ($\alpha\to 0$),
the transmission is equivalent to the case of single-photon
scattering.Sánchez-Burillo  _et al._ (2016) The superscript $t^{(1)}$ denotes
that we are considering the single molecule scenario. A single molecule can be
treated as a point-like object, leading to the continuity of the photon
wavepacket at the interaction point with the molecule and the guide, as
expressed by:

$r^{(1)}(\omega)=t^{(1)}(\omega)-1.$ (36)

For the calculation of transmission (Eq. (35)), we refer to Eq. (30). The two-
level system (TLS) dynamics are described by the master equation:

$\displaystyle\frac{d\varrho}{dt}=$ $\displaystyle-i[{\cal
H}_{s},\varrho]-i\alpha\sqrt{\frac{d\omega}{dk}}g(\omega_{0})\cos(\omega_{0}t)[\sigma_{x},\varrho]$
(37)
$\displaystyle+2\gamma(\Omega)\left(\sigma^{-}\varrho\sigma^{+}-\frac{1}{2}\\{\sigma^{+}\sigma^{-},\varrho\\}\right)$
$\displaystyle+2\gamma(\Omega)e^{-\beta\Omega}\left(\sigma^{-}\varrho\sigma^{+}-\frac{1}{2}\\{\sigma^{+}\sigma^{-},\varrho\\}\right)$
$\displaystyle+\gamma_{\phi}\left(\sigma^{+}\sigma^{-}\varrho\sigma^{+}\sigma^{-}-\frac{1}{2}\\{\sigma^{+}\sigma^{-},\varrho\\}\right)\,.$

Here,

${\cal H}_{s}=\Omega\sigma^{+}\sigma^{-},\quad\Omega\equiv g\mu_{B}B$ (38)

and $\gamma$ represents the temperature-dependent rates for the spins,
including both the coupling to the line and any intrinsic loss mechanisms
within the molecules:

$\gamma(\Omega)=\left(\xi+2\pi\lambda_{\Omega}^{2}\right)(n_{\Omega}+1),$ (39)

with $n_{\Omega}=1/(e^{\beta\Omega}+1)$ being the Bose occupation number, and
$\xi$ accounting for intrinsic loss mechanisms. Additionally, $\gamma_{\phi}$
quantifies the pure dephasing terms. The factor of $2$ in front of the
summation is because the master equation is coupled to both left
($l_{\omega}$) and right ($r_{\omega}$) modes. Solving the master equation
within the linear response theory (LRT), we derive:

$t^{(1)}(\omega)=1-\frac{2\pi\lambda_{\omega}^{2}\langle\sigma_{z}\rangle_{\beta}}{\Gamma+i(\Omega-\omega)},$
(40)

where

$\displaystyle\langle\sigma_{z}\rangle$
$\displaystyle=\tanh\left(\frac{\beta\omega}{2}\right),$ (41)
$\displaystyle\Gamma$ $\displaystyle=\lambda(\Omega)+\gamma_{\phi}.$ (42)

Figure S16: Graphical representation of the transfer matrix, $T$.

It is essential to highlight that the transfer matrix, which describes the
relationship between the photon wavepacket on both sides of the molecule
(refer to Fig. S16), can be expressed in terms of the transmission and
reflection coefficients as follows:

$T=\begin{pmatrix}t-r^{2}/t&r/t\\\ -r/t&1/t\end{pmatrix}.$ (43)

Using Eq. (36), the transmission for a single molecule can be succinctly
expressed as:

$T^{(1)}=1+\frac{r^{(1)}}{t^{(1)}}(\sigma_{x}+i\sigma_{y}).$ (44)

To prevent any misunderstanding, it is worth emphasizing that the Pauli
matrices mentioned above do not pertain to those representing $1/2$-spins. In
this context, they serve merely as a convenient method for representing a
$2\times 2$ matrix.

###  VI.2 Transmission in the paramagnetic phase

In the paramagnetic phase, the molecules act as independent emitters.
Consequently, the total transfer matrix is determined by the product of
individual matrices,

$T=\prod
T_{j}^{(1)}=1+\sum_{j}\frac{r^{(1)}_{j}}{t^{(1)}_{j}}(\sigma_{x}+i\sigma_{y}),$
(45)

leveraging the property of the Pauli matrices,

$(\sigma_{x}+i\sigma_{y})^{2}=0.$ (46)

It should be noted that the total transmission is found as [cf. Eq. (43)],

$t_{\text{total}}=\frac{1}{1+\sum_{j}\frac{r^{(1)}_{j}}{t^{(1)}_{j}}}.$ (47)

By employing Eqs. (40) and (36), and recognizing that
$\lambda_{\omega}\ll\Gamma$ in Eq. (42), we find that

$\frac{r^{(1)}_{j}}{t^{(1)}_{j}}\approx\frac{2\pi\lambda_{\omega}^{2}\langle\sigma_{z}\rangle}{-\Gamma+i\delta}\Longrightarrow\sum_{j}\frac{r^{(1)}_{j}}{t^{(1)}_{j}}\approx\frac{N2\pi\lambda_{\omega}^{2}\langle\sigma_{z}\rangle}{-\Gamma+i\delta},$
(48)

where

$\delta=\Omega-\omega.$ (49)

Thus, under these assumptions, the transfer matrix across the ensemble behaves
as though scattering were occurring through a single object, but with enhanced
coupling $\sum_{n}g_{j}^{2}$, i.e., superradiance. In the main text, we define
the collective coupling as

$G=N2\pi\lambda_{\omega}^{2}\langle\sigma_{z}\rangle$ (50)

leading to

$S_{21}=\left|t_{\rm
total}\right|=\left|\frac{1}{1-\frac{G}{-\Gamma+i\delta}}\right|,$ (51)

which is Eq. (1) in the main text.

###  VI.3 Transmission below $T_{N}$

When intrinsic spin interactions become significant, DPPH molecules tend to
organize into AFM 1D chains. Consequently, microwave photons no longer
interact with individual emitters but with spin waves ($k=0$) that behave as
bosons. We know that transmission through a boson is given by [cf with (40)],

$t^{(1)}(\omega)=1-\frac{2\pi\lambda_{\omega}^{2}}{\Gamma_{\rm
SW}+i(\Omega-\omega)}\;.$ (52)

It is crucial to realize that, now, $\Gamma_{\rm SW}$ is not temperature-
dependent, meaning that transmission does not depend on temperature.
Furthermore, due to anisotropy, each chain possesses a distinct angle $\psi$,
which influences its resonance frequency as we have calculated and illustrated
in Figure S15. Thus, the sum over $j$ in (51) is across chains, each with a
different frequency $\Omega(\psi)$. We can replace the sum with an integral,
resulting in

$S_{21}=\left|t_{\rm total}\right|=\left|\frac{1}{1+\int
d\psi\sin\psi\frac{\lambda_{\Omega}^{2}N_{\rm
eff}(T_{N})}{-\Gamma_{SW}+i\delta_{\psi}}}\right|\ ,$ (53)

with $\Gamma_{SW}\cong\Gamma(T_{N})$ and $\delta_{\psi}=\Omega(\psi)-\omega$,
which is Eq. (4) in the main text.

## References

  * Williams (1966) D. E. Williams, J. Am. Chem. Soc. 88, 5665 (1966).
  * Žilić  _et al._ (2010) D. Žilić, D. Pajić, M. Jurić, K. Molčanov, B. Rakvin, P. Planinić, and K. Zadro, J. Magn. Reson. 207, 34 (2010).
  * Yordanov (1996) N. D. Yordanov, Appl. Magn. Reson. 10, 339 (1996).
  * Kent  _et al._ (1994) A. D. Kent, S. von Molnár, S. Gider, and D. D. Awschalom, J. Appl. Phys. 76, 6656 (1994).
  * Ohya-Nishiguchi (1979) H. Ohya-Nishiguchi, Bull. Chem. Soc. Jpn. 52, 3480 (1979).
  * Fujito (1981) T. Fujito, Bull. Chem. Soc. Jpn. 54, 3110 (1981).
  * Stewart (1983) G. R. Stewart, Rev. Sci. Instrum. 54, 1 (1983).
  * Blöte (1975) H. Blöte, Physica B+C 79, 427 (1975).
  * Vleck (1948) J. H. V. Vleck, Phys. Rev. 74, 1168 (1948).
  * Anderson and Weiss (1953) P. W. Anderson and P. R. Weiss, Rev. Mod. Phys. 25, 269 (1953).
  * Gimeno  _et al._ (2020) I. Gimeno, W. Kersten, M. C. Pallarés, P. Hermosilla, M. J. Martínez-Pérez, M. D. Jenkins, A. Angerer, C. Sánchez-Azqueta, D. Zueco, J. Majer, A. Lostao, and F. Luis, ACS Nano 14, 8707 (2020).
  * Hümmer  _et al._ (2013) T. Hümmer, F. J. García-Vidal, L. Martín-Moreno, and D. Zueco, Phys. Rev. B 87, 115419 (2013).
  * Ghirri  _et al._ (2016) A. Ghirri, C. Bonizzoni, F. Troiani, N. Buccheri, L. Beverina, A. Cassinese, and M. Affronte, Phys. Rev. A 93, 063855 (2016).
  * Mergenthaler  _et al._ (2017) M. Mergenthaler, J. Liu, J. J. Le Roy, N. Ares, A. L. Thompson, L. Bogani, F. Luis, S. J. Blundell, T. Lancaster, A. Ardavan, G. A. D. Briggs, P. J. Leek, and E. A. Laird, Phys. Rev. Lett. 119, 147701 (2017).
  * Garcia-Palacios  _et al._ (2009) J. Garcia-Palacios, J. Gong, and F. Luis, J. Phys. Condens. Matter 21, 456006 (2009).
  * Vonsovskii (1966) S. V. Vonsovskii,  _Ferromagnetic resonance_ (Pergamon Press, 1966).
  * Ross (2013) M. P. Ross,  _Spin dynamics in an antiferromagnet_ , Diploma thesis, Technische Universität München (2013).
  * Gardiner and Collett (1985) C. W. Gardiner and M. J. Collett, Phys. Rev. A 31, 3761 (1985).
  * Sánchez-Burillo  _et al._ (2016) E. Sánchez-Burillo, L. Martín-Moreno, J. J. García-Ripoll, and D. Zueco, Phys. Rev. A 94, 053814 (2016).