\section{Entanglement capacity}

The entanglement capacity $C_{E}$ \cite{Yao2010,Schliemann2011,Boer2019,ENT:Nandy2021EntanglementCapacity} that we define in the main text is a measure of how susceptible the entanglement stored in the system is upon its extraction. It can be interpreted as the quantum information theoretic counterpart of the thermodynamic heat capacity. It is obtained starting from the entanglement Hamiltonian $\mathcal{K}_A\equiv-\ln\rho_A$, with $\rho_A$ being the density matrix after tracing-out rest of the system in a given bipartition $A$.  The entanglement capacity is defined as:
\begin{equation}\label{EqA:EntCap}
C_E(\rho_A)=\text{Tr}(\rho_A\mathcal{K}_A^2)-\text{Tr}(\rho_A\mathcal{K}_A)^2=\text{Tr}(\rho_{A}(\ln\rho_{A})^2)-(\text{Tr}(\rho_{A}\ln\rho_{A}))^2.
\end{equation}

In particular the entanglement capacity can be related to the Renyi entropies $S^{(m)}_A=\frac1{1-m}\ln\text{Tr}\rho_A^m$ by
\begin{equation}\label{EqA:EntCapRenyi}
C_E(\rho_A)=\lim_{m\rightarrow1}m^2\partial_m^2((1-m)S^{(m)}_A).
\end{equation}

\section{Details about Numerical Simulations}

In order to study the criticality of the light-matter system, we perform numerical simulations using exact diagonalization (ED)~\cite{Sandvik2010} and density matrix renormalization group (DMRG)~\cite{white_prl_1992,white_prb_1993, schollwock_aop_2011, Orus_aop_2014} calculations. 

Within the ED simulations, the Hamiltonian of the system is written in a sparse matrix form and a Lanczos algorithm is used to find the ground state and the first few excited states. We only consider the system at half-filling so that the Hilbert space dimension for the matter part is $L!/(L/2)!^2$, corresponding to $\sim10^5$ for $L=20$. In principle, the photon Hilbert space is infinite dimensional, but in practice it can be truncated up to a maximum photon number $n^{\max}_{\textrm{ph}}$. The total dimension of the system Hilbert space makes ED calculations viable only for small system sizes ($L\lesssim20$) with limited maximum number of photons $n^{\max}_{\textrm{ph}}\sim 5$. It is to be noted that the average photon occupancy $n_{\textrm{ph}} = \braket{\hat a^{\dagger} \hat a}$ increases linearly with system-size in the symmetry-broken charge-density wave (CDW) phase (see Fig. 2(c) of the main text); thus due to the small value of $L$, the photon population is not too large and typically saturates to $n_{\textrm{ph}}\lesssim 3$. Despite the large computation cost, ED calculations have the advantage of giving access to the full quantum state, allowing the study of any degeneracy of the ground state.

\begin{figure}[htb]
    \centering
    \includegraphics[width=0.5\linewidth]{Omega_1_conv.pdf}
    \caption{(a) The average photon occupancy $n_{\textrm{ph}}$ and the electric field $E = |\braket{\hat a + \hat a^{\dagger}}|$ as functions of the maximum photon cut-off $n^{\max}_{\textrm{ph}}$ in the symmetry broken CDW phase ($g_{\textrm{eff}}/t=1.5$) for $\omega/t=1$ and system-size $L=100$. Both $n_{\textrm{ph}}$ and $E$ saturates to converged values for $n^{\max}_{\textrm{ph}} \gtrsim 70$.
    (b) The ground state energy difference $E_0(\chi) - E_0(\chi=600)$ with respect to the ground state energy obtained with MPS bond dimension $\chi=600$ and the total DMRG truncation error with varying MPS bond dimesion $\chi$. Here, we set $\omega/t=1$ and $L=100$, and consider $g_{\textrm{eff}}/t=0.962$ where the photon entanglement $S_P$ is highest.
    }
    \label{fig:conv}
\end{figure}


The DMRG calculations are based on a matrix-product states (MPS)~\cite{schollwock_aop_2011, Orus_aop_2014} ansatz for a finite system size $L$ with open boundary condition and for truncated photon space. Compared to ED simulations, DMRG calculations allow us to access larger system sizes ($L \sim 100$) with a higher cut-off in the photonic Hilbert-space ($n^{\max}_{\textrm{ph}} \sim 100$), enough to ensure convergence with respect to the truncated photon space (see Fig.~\ref{fig:conv}(a)). Thus, our MPS ansatz consists of $L+1$ sites where the first MPS site corresponds to the photon with local dimension $d_{\textrm{ph}} = n^{\max}_{\textrm{ph}}+1 = 101$, while the other sites are for the fermions with local dimension of $2$. Moreover, to keep the fermion number fixed, we utilize $U(1)$ symmetry-preserving tensors \cite{Singh_PRA_2010,Singh_PRB_2011} at the fermionic MPS sites. Such an MPS setup has been proven to be very successful in similar hybrid systems~\cite{Halati20a, Halati20b, Halati21}. The maximal bond dimension for the MPS representation is set to $\chi=600$, which is sufficient to keep the total discarded weight of the singular values below $10^{-10}$ ($\sim 10^{-12}$ per MPS bond) and ensure convergence for simulations with system sizes up to $L=100$ (see Fig.~\ref{fig:conv}(b)).

\section{Results for $\omega/t=3$ and $\omega/t=10$}

In the main text, we have presented our analysis for $\omega/t=1$ only. To supplement those results and for the purpose of completeness,
in this section, we present the same analysis for $\omega/t=3$ (Fig. \ref{figS:collapse_Om3}) and $\omega/t=10$ (Fig. \ref{figS:collapse_Om10}).

We find that the results for $\omega/t=3$ are qualitatively similar to those obtained for $\omega/t=1$: in particular, the collapse of the finite size scaling analysis is very good and the critical exponents for entanglement and order parameter are in excellent agreement. On the other hand the quality of the collapse is worse for $\omega/t=10$ as seen by the larger errors and by the fact that the differences in the critical exponent $\nu$ and in the critical coupling $g_{\textrm{eff}}^c/t$  extracted from entanglement and order parameter are also larger.

\begin{figure*}[!t]
    \centering
    \includegraphics[width=\textwidth]{Omega_3.pdf}
    \caption{Phase transition analysis at $\omega/t=3$ for the data obtained from DMRG simulations. (a) The photon entanglement entropy $S_P$ and the order parameter $O_D$ as functions of the effective coupling $g_{\textrm{eff}}$. (b) Collapse of the curves using the scaling hypothesis near the transition point. As in the main text, we consider system sizes in the range $L \in [48, 100]$ for the scaling analysis, and plot the data only for $L=48, 64, 80, 96$ for clarity. The values of the critical coupling and of the critical exponents for the two transitions are reported in the box; we observe that the values of $g^c_{\textrm{eff}}$ and $\nu$ agree within error. (c) Entanglement capacity $C_E(\rho_P)$ (solid lines) and photon entropy $S_P$ (dashed). (d) Scaling of the locations $g_{\textrm{eff}}^{\max}(L)$ of the peak values with respect to the system size $L$ for both $S_P$ and $C_E(\rho_P)$. For both quantities, $g_{\textrm{eff}}^{\max}$ approaches the critical point $g_{\textrm{eff}}^c$ in the thermodynamic limit with a power law behavior.}
    \label{figS:collapse_Om3}
\end{figure*}

We postulate that this difference is model dependent: one possible explanation is due to the higher photon energy which makes it difficult to excite and entangle with the electrons. We also remark that in the absence of coupling, the spectrum of the electrons has a bandwidth of $4t$, so that no direct excitations can be mediated by the photon when $\omega \gtrsim 4t$. Thus, we may expect qualitative differences when $\omega/t \gtrsim 4$.

\begin{figure*}[!t]
    \centering
    \includegraphics[width=\textwidth]{Omega_10.pdf}
    \caption{Phase transition analysis at $\omega/t=10$ for the data obtained from DMRG simulations. (a) The photon entanglement entropy $S_P$ and the order parameter $O_D$ as functions of the effective coupling $g_{\textrm{eff}}$. (b) Collapse of the curves using the scaling hypothesis near the transition point. Again we consider system sizes in the range $L \in [48, 100]$ for the scaling analysis, and plot the data only for $L=48, 64, 80, 96$. The values of $g^c_{\textrm{eff}}$ and $\nu$ still agree within uncertainty, but the difference is now larger compared to the $\omega/t=1$ and $\omega/t=3$ cases; the uncertainties on the fit parameters are also larger, signalling a worse quality of the collapses. (c) Entanglement capacity $C_E(\rho_P)$ (solid lines) and photon entropy $S_P$ (dashed). (d) Scaling of the locations $g_{\textrm{eff}}^{\max}(L)$ of the peak values with respect to the system size $L$ for both $S_P$ and $C_E(\rho_P)$. The asymptotic values of $g_{\textrm{eff}}^{\max}$ are again compatible with $g_{\textrm{eff}}^c$.}
    \label{figS:collapse_Om10}
\end{figure*}

We also report in Fig. \ref{figS:electric field} the scaling analysis and the collapsed curves for the rescaled electric field $E/\sqrt{L}$ for $\omega/t=1$. The uncertainty on the critical parameters are very small and the collapses are excellent, signalling the high quality of the fit. We remark that the electric field exhibits the same critical effective coupling and the same value of $\nu$ as the entanglement entropy of the photon and the order parameter, meaning it can also be used to extract the critical exponents experimentally.

\begin{figure*}[!t]
    \centering
    \includegraphics[width=0.6\textwidth]{Omega_1_current.pdf}
    \caption{Critical behavior of the rescaled electric field $|\langle\hat a+\hat a^{\dagger}\rangle|/\sqrt L$ at $\omega/t=1$ for data obtained from DMRG simulations. (a) Plot of the rescaled electric field as functios of the effective coupling $g_{\textrm{eff}}$ for selected system sizes. (b) Collapse of the curves  obtained from a finite size scaling analysis of system-sizes in the range $L \in [48, 100]$ for the scaling analysis, and are plotted only for $L=48, 64, 80, 96$ for clarity. The values of the critical coupling and of $\nu$ coincide with those obtained for the entanglement entropy and the order parameter.}
    \label{figS:electric field}
\end{figure*}


\section{Choice of the symmetry breaking term}

In the main text, we have chosen a symmetry breaking term of the form $g_2(\hat c^{\dagger}_1\hat c_1-\hat c^{\dagger}_L\hat c_L)$, which selects only one of the possible charge density wave states. There is some degree of arbitrariness in this choice: for example one could choose a term that breaks the degeneracy in the electric field, such as $g_2(\hat a+\hat a^{\dagger})$. Indeed, we have performed calculations employing such symmetry breaking term and found results qualitatively similar to those presented in the rest of the work. However, the main drawback of such choice is that the expectation value of the electric field $\hat E=\hat a+\hat a^{\dagger}$ does not depend on the system size $L$ in the disordered phase and scales like $\sqrt{L}$ in the ordered phase. This implies that the symmetry breaking term scales differently with $L$ depending on the value of $g$, making a finite size scaling analysis imprecise. Indeed, we have compared the collapses obtained from the scaling fit and observed that the ones obtained from $g_2(\hat c^{\dagger}_1\hat c_1-\hat c^{\dagger}_L\hat c_L)$ are significantly better.


\end{document}