Microcavity phonon polaritons from the weak to the ultrastrong phonon–photon coupling regime

Strong coupling between molecular vibrations and microcavity modes has been demonstrated to modify physical and chemical properties of the molecular material. Here, we study the less explored coupling between lattice vibrations (phonons) and microcavity modes. Embedding thin layers of hexagonal boron nitride (hBN) into classical microcavities, we demonstrate the evolution from weak to ultrastrong phonon-photon coupling when the hBN thickness is increased from a few nanometers to a fully filled cavity. Remarkably, strong coupling is achieved for hBN layers as thin as 10 nm. Further, the ultrastrong coupling in fully filled cavities yields a polariton dispersion matching that of phonon polaritons in bulk hBN, highlighting that the maximum light-matter coupling in microcavities is limited to the coupling strength between photons and the bulk material. Tunable cavity phonon polaritons could become a versatile platform for studying how the coupling strength between photons and phonons may modify the properties of polar crystals.

1 Description of the cavities and permittivities of the different materials In this Supplementary Note, we describe in more detail the hexagonal boron nitride (hBN) microcavities that we analyze in this work, including all the different thicknesses. We use Fabry-Pérot cavities containing hBN layers of variable thickness. A schematic diagram is shown in Supplementary   Fig. 1. The cavities are formed by planar layers, and all the interfaces are perpendicular to the z axis, as shown by the coordinates axes included in the scheme. In all the calculations the incident medium (i.e. the medium from which the system is illuminated) is vacuum, the substrate is CaF 2 and the mirrors are 20 nm-thick gold layers. The inside of the cavity, i.e. without considering the mirrors, extends from z = 0 to z = L cav (L cav is the total thickness), and it contains a layer of thickness L hBN placed between z = L 1 and z = L 2 . In all the experiments and in many of the transfer-matrix simulations, the material of this layer is hBN, and we refer to this particular case as a hBN-filled cavity. In other simulations, we replace hBN by a material with constant permittivity ε hBN,∞ , which corresponds to the high-frequency permittivity of hBN (see below). We refer to these cavities as bare cavities, because the contribution of the hBN phonons is eliminated. The rest of the cavity (the layers between z = 0 and z = L 1 and between z = L 2 and z = L cav ) is filled by MoS 2 , which does not show any resonant feature in the analyzed range of frequencies.
In the transfer-matrix simulations, we have modelled the materials of the system as follows. For the incident medium we use the vacuum permittivity, and the relative permittivity of the substrate material CaF 2 is ε CaF 2 = 1.882. The gold mirrors are described by a Drude permittivity that fits the low-frequency experimental data from Supplementary Ref. [1]: with plasma frequency ω p,Au = 73114.15 cm −1 and damping frequency γ Au = 571.04 cm −1 [2].
We model hBN with a diagonal permittivity tensor ↔ ε hBN (ω), in order to take into account the anisotropy of this material. hBN is characterized by a layered atomic structure, so that the phonons that are polarized along the plane of the atomic layers (in-plane phonons) have different frequencies compared to those polarized in the normal direction, i.e. parallel to the anisotropic axis ( , out-ofplane phonons). In our configuration, the plane of the atomic layers is parallel to all the Fabry-Pérot interfaces (x − y plane). We thus have ↔ ε hBN (ω) = diag ε hBN (ω), ε hBN (ω), ε hBN, (ω) , where we omit the second subscript in the x and y components for brevity (notice that the x and y components are identical). In particular, the in-plane permittivity tensor component ε hBN (ω) follows a Lorentzian function [3]: where ω TO is the corresponding transverse optical (TO) phonon frequency, ω LO is the longitudinal optical (LO) phonon frequency, γ is the damping frequency and ε hBN,∞ is the high-frequency permittivity. For our analysis, the exact value of ω TO determines the detuning between the cavity mode and the phonon, which is a crucial parameter in strong coupling. Therefore, to obtain the value of this parameter in our samples, we measured the reflectivity spectra of a 100 nm-thick hBN layer before the fabrication of the cavities. A clear resonance was observed at ω TO = 1364 cm −1 , which is the value that we consider in the simulations. This frequency is very close to the value ω TO = 1360 cm Ref. [3]).
In order to obtain the permittivity of MoS 2 , we have fabricated three microcavities fully filled by MoS 2 of different nominal thickness: L cav = 660 nm, 740 nm and 960 nm. The geometry is the same as in Supplementary Fig. 1 but without the hBN layer. We measure the reflectivity spectra and obtain the frequencies of the Fabry-Pérot modes for each cavity up to 8000 cm −1 . We then perform transfer-matrix simulations of the cavities and find the value of the permittivity ε MoS 2 required to reproduce the position of the experimental dips. We show in Supplementary Fig. 2 where ω is in units of cm −1 .
Last, in some simulations in the end of the main text, we replace hBN with PMMA and a weak oscillator. For PMMA, we use a permittivity that also follows a Lorentzian form [4]: 2 Characterization of the polaritonic frequencies

Transfer-matrix simulations
We first summarize briefly how to calculate the reflectivity spectra of hBN microcavities using the transfer-matrix formalism [5]. This method allows to obtain the total reflection r s(p) total and transmission t s(p) total coefficients of a s(p)-polarized planewave, for any layered structure, in terms of the thickness d i and relative permittivity ε i of each layer i and the Fresnel coefficients r s(p) i,i+1 and t s(p) i,i+1 of each interface between layers i and i + 1. In total, the system contains N lay finite-sized layers, and we label the incident medium as i = 0 and the substrate as i = N lay + 1. For the calculations in the main text, we consider that the light is incident along the z direction perpendicular to the planar interfaces and all materials are considered isotropic. However, for the calculations in Supplementary Notes 4 and 5, we analyze the effects of the anisotropy of hBN and, thus, here we discuss the complete formalism for anisotropic materials. To take the anisotropy of hBN into account, we consider that each layer i has a permittivity tensor of the general diagonal form ↔ ε i = diag(ε i,x , ε i,y , ε i,z ) (ε hBN,x = ε hBN,y for hBN in our experiments). Light is incident in the x − z plane. The Fresnel coefficients are given by [5] where ξ = k x /k 0 is the parallel component of the wavevector k x normalized with respect to the wavevector in vacuum k 0 .
The total transfer matrix T s(p) relates the amplitudes of the electromagnetic field in the incident medium (i = 0) and the substrate (i = N lay +1) for s(p)-polarized light: A A s i,± and A p i,± correspond to the amplitude of the electric and magnetic field, respectively (the amplitude of the field parallel to the interface for each polarization). The second subindex indicates the direction of propagation. + corresponds to the direction of the incoming planewave (and thus also to the direction of the transmitted light) and, similarly, the subindex -corresponds to the direction of the reflected light. In our system, A s(p) N lay +1,− = 0, because in the substrate the electromagnetic field only propagates in the direction of the transmission. T s(p) is given by the expression with The matrix M s(p) i,i+1 relates the amplitudes of the electromagnetic field propagating both in the + and − directions, at both sides of the interface between layers i and i + 1. is obtained, these matrix elements are used to calculate the total transmission and reflection coefficients of the system as t and r s(p) . Last, we obtain the reflectivity spectra R s(p) , which is defined as the ratio between the intensity of the reflected and incident light. It is obtained as R s(p) = |r total corresponds to the ratio of the amplitude of the electromagnetic field). Once r s(p) total and t s(p) total are known, this approach can also be used to calculate the electromagnetic fields in any point of space. We stress that, as discussed above, in every transfer-matrix simulation shown in this work except in Supplementary Notes 4 and 5, we consider light incident in the z direction (ξ = 0). Under this condition, the z component of the permittivity does not affect the results and, taking into account that ε hBN,x = ε hBN,y for hBN, we recover the usual equations for isotropic materials. The anisotropy needs to be considered in Supplementary Note 4 for focused illumination and in Supplementary Note 5 for oblique-incident illumination.

Frequencies of the modes
The coupling strength g is determined from the position of the eigenfrequencies of the system. In general, however, the reflectivity dips (or transmission peaks) that are found in experimental optical spectra do not need to coincide with the eigenfrequencies, particularly for small coupling strengths. To characterize these complex eigenfrequencies, we also use the transfer-matrix formalism. From the total transfer matrix T s(p) , we obtain the analytical expression for the reflection and transmission coefficients r s(p) total and t s(p) total of the total system. Then, we solve numerically the complex frequencies for which the denominator of these coefficients vanishes (i.e. the poles), by extending all the analytical expressions for the permittivities (see Supplementary Note 1) to complex values of ω. a The eigenfrequencies of the bare-cavity modes and of the polaritonic modes correspond to the complex poles of the reflection and transmission coefficients of the bare and hBN-filled cavities, respectively. The real part of these eigenfrequencies is the mode frequency and the imaginary part is equal to (minus) half the losses of the mode.
We compare in Supplementary Fig. 4 and in Figs. 3 and 5 of the main text the real part of these eigenfrequencies (indicated by dashed lines) with the frequencies of the reflectivity dips for the three cavities explored in more detail in this work (100 nm-thick hBN layer, 10 nm-thick hBN layer and the fully filled cavity, respectively). The agreement between these two frequencies is very good, and thus in these systems the reflectivity dips are a good indication of the bare-cavity and polaritonic modes.
However, the calculation of the eigenfrequencies is important for determining the coupling strength for very thin layers in Fig. 4 of the main text.

Coupled harmonic oscillator model
In this subsection, we explain the procedure that we follow in the main text to compute the coupling strength g. This method is based on modelling the interaction between the material and the cavity mode by means of two coupled harmonic oscillators. This discussion is complemented in Supplementary Note 6 with an alternative method to calculate g, where we derive an analytical expression based on the microscopic interaction of the phononic material with the electric field of the cavity mode.
In this model, one oscillator (x (j) cav ) is associated to the electromagnetic field of the mode of order j and frequency ω (j) cav of the bare cavity, a cavity where we eliminate the phonon contribution and describe the permittivity of hBN as ε hBN,∞ . The second (x TO ) represents the phonon at frequency ω TO (which is modelled in Supplementary Note 6 explicitly as the collective excitation of the dipoles associated with the atomic vibrations at each unit cell within hBN). The coupling strength between these two oscillators is g, and the equations of motion are [6][7][8][9]: κ and γ are the losses of the cavity and the phonon, respectively, and the dot · denotes the time derivative. We note that we use the derivative of x TO and x (j) cav in the terms describing the coupling (last term in the left hand side of both equations). This choice is necessary to reproduce the results obtained in the ultrastrong coupling regime with transfer-matrix simulations and experimentally [9]. a As we use simple Lorentzian and Drude functions, the expression of the permittivities can be directly evaluated at complex frequencies. In our system, the losses of the cavity vary between κ ≈ 30 cm −1 , for the cavity fully filled with hBN, and κ ≈ 60 cm −1 , for cavities with very thin layers of hBN. These values yield quality factors Q = ω (j) cav κ between Q ≈ 22 and Q ≈ 45 for all cavities. On the other hand, the experimental measurement of this value κ is challenging. However, by comparing the widths of the measured reflectivity polaritonic dips with the simulations, we observe that the experimental widths are only slightly larger, and we estimate that the difference between the experimental and theoretical value of κ is at most 10% (notice that the sum of the widths of the two polaritonic dips should be close to κ + γ [10]).
Supplementary equation (9) can be alternatively written in frequency domain: The polaritonic frequencies ω (j) ± correspond to the values that cancel the determinant of this system of equations, i.e., which can be solved numerically for any value of the coupling strength. Further, approximate analytical solutions can be obtained for two typical situations. For small values of g, the polaritonic frequencies are similar to the bare frequencies: ω cav , ω TO . Thus, we can make the secular approximation (ω (and the corresponding approximation for ω TO ) [11]. Then, Supplementary Eq. (10) becomes and Supplementary Eq. (11) becomes quadratic in ω On the other hand, Supplementary Eq. (11) can also be solved analytically for the limiting case where the losses κ and γ are negligible compared to g. Under this condition, the polaritonic frequencies are To calculate the coupling strength g, we first extract the polaritonic frequencies ω H Hop [12]. In this alternative approach, light and matter excitations are also represented by harmonic oscillators which interact with coupling strength g, but the behaviour of the system is described by the Hamiltonian whereâ andâ † are the annihilation and creation operator of the j th cavity mode, respectively, whileb L cav (nm) 500 1000 andb † are the corresponding operators for the TO phonon. Apart from the interaction term proportional to g, this Hamiltonian also includes the diamagnetic term, which is proportional to g 2 [13,14].
Due to the diamagnetic and counter-rotating terms ofĤ Hop , the energy of the ground state is shifted from zero. The shift of the ground state is not described by the classical harmonic oscillator model, but in this study we are interested in the transition energies between the excited states and the ground state. The result obtained for these transition energies using the Hopfield Hamiltonian is equal to the polaritonic energies obtained from Supplementary Eq. (14). Thus, the classical harmonic oscillator model is able to describe the transition energies also in the ultrastrong coupling regime [9].
Last, we note that the harmonic oscillator model presented here treats each cavity mode of order j independently, i.e. that different cavity modes do not interact with each other and thus the coupling between hBN and each cavity mode is described by a separate Supplementary Eq. (9). For the cavities that we analyze in this work, which are designed so that the TO phonon interacts mostly with the first cavity mode, the results in Supplementary Fig. 3 show that this approach is appropriate. However, we have found that for thicker cavities where higher-order modes are tuned to the TO phonon, the harmonic oscillator does not always describe well the position of the polaritonic frequencies as obtained from the transfer-matrix simulations. We believe that solving this discrepancy requires to include in the harmonic oscillator model the possibility of coupling different modes with each other, for instance following an approach based on quasinormal modes [15,16] (a similar reason may be behind the small discrepancy between the two values of the coupling strength that are obtained in Supplementary Fig.   8 for the 3 rd cavity mode, as discussed in Supplementary Note 6). However, the results of the cavities analyzed in this work would not be significantly affected by future analysis of this type. 3 Reflectivity spectra of the cavity embedding a 100 nm-thick hBN layer In the main text, we show the reflectivity spectra of a cavity embedding a 10 nm-thick hBN layer (Fig. 3) and a cavity fully filled by hBN (Fig. 5), as a function of the total cavity thickness L cav .
To complement these results, we show in Supplementary Fig. 4a the reflectivity spectra of the cavity embedding a 100 nm-thick hBN layer. For reference, in Supplementary Fig. 4b we show a calculation of a bare cavity where the 100 nm layer is formed by a material with constant permittivity ε hBN,∞ .
The dashed lines in Supplementary Fig. 4a indicate the polaritonic frequencies ω  Supplementary Fig. 4b correspond to the bare cavity mode frequencies ω (j) cav . As expected, we observe that the anticrossing between the polaritonic frequencies for the 100 nm-thick hBN layer is considerably smaller than for the fully filled cavity, but much larger than for the cavity with a 10 nm layer.

Reflectivity of the system under focused illumination
In the transfer-matrix calculations of the main text, we assume that the propagation direction of the incident light is normal to the surface of the mirrors and the substrate (z direction). However, the experiments have been performed using a microscope with focused illumination. In this Supplementary Note, we discuss the effect of this focusing on the reflectivity spectra of the system. We show that this effect is small and that the main difference with respect to the spectra under normal incident light is the possibility to observe the coupling with the out-of-plane phonons of hBN, instead of just with the in-plane phonons (see Supplementary Note 1 for a discussion of these two types of phonons).
The total reflectivity R foc for focused illumination is obtained by decomposing the incident light as an integral over planewaves incident at different angles. We then obtain (for more details, see Chapter 3.9 in Supplementary Ref. [17]): where I in and I ref correspond to the intensity of the incident and reflected focused beam, respectively, total is the ξ-dependent reflection coefficient of the full system for a s(p)-polarized planewave, which is discussed thoroughly in Supplementary Note 2.1. The upper limit of the integral is given by the numerical aperture of the microscope used to focus light, NA = sin θ max . In our case, we have NA = 0.4, and thus, θ max ≈ 23.5 • . Furthermore, we set θ min = 10 • for the lower bound of the integral, because the microscope used in the experiments obstructs the propagation of the central part of the light incoming to the focusing lens, which eliminates the contribution from the small-angle components.
We have verified that, for cavities embedding a 10 nm layer of hBN, the results of the transfermatrix simulations using focused illumination are nearly identical compared to the spectra for normal incidence. Crucially, light travels at considerably smaller angles inside this cavity than in free space, due to the high permittivity of MoS 2 , ε MoS2 , which fills most of the cavity for such a thin hBN layer.
In particular, the angles inside MoS 2 range from 2.6 • to 5.9 • , and thus the difference with the case of normal incidence is very small. We next consider a fully filled hBN cavity of thickness L cav = 1665 nm (whose first cavity mode is resonant with the TO phonon frequency) and show in Supplementary Fig. 5(a-c) the reflectivity spectra under focused (red dashed line) and normal-incidence illumination (black solid line), both calculated by transfer-matrix simulations. The effect of using focused light becomes larger for the cavity fully filled by hBN, as the high-frequency permittivity of hBN ε hBN,∞ is smaller than ε MoS2 , and thus the angle of light propagation inside the cavity can be larger. A small but appreciable difference between normal incident light and focused illumination is observed for high wavenumbers, as shown in panel a. The reflectivity dips obtained under focused illumination are displaced towards larger wavenumbers compared to the normal-incidence spectra, because the frequencies of the Fabry-Pérot cavity modes depend on the angle of propagation. Furthermore, the widths of the dips are larger for focused illumination because they are the result of the sum of different contributions, each corresponding to a different angle and thus resonant at a slightly different wavenumber (Supplementary Eq. (16)). However, these differences remain small, and the agreement between the two spectra improves further for frequencies close to ω TO (see zoom in panel b), which is the main region of interest in this work. Thus, the use of normal incidence in our calculations of the other sections is justified.
Possibly the most interesting feature of the results obtained under focused illumination is observed in Supplementary Fig. 5c, which shows the reflectivity spectra at low wavenumbers (near the frequencies of the out-of-plane phonons). An additional reflectivity dip for focused illumination (red dashed line) appears compared to the normal-incidence spectra (black solid line), which is due to the anisotropy of hBN. For normal incidence, the illumination only couples with the in-plane transverse optical phonon, which is polarized in the parallel direction to all planar interfaces. On the other hand, for focused illumination it becomes possible to excite the phonons that are polarized in the out-of-plane direction and that are found at significantly lower frequencies (ω TO, = 746 cm −1 and ω LO, = 819 cm −1 , the Reststrahlen band of the out-of-plane phonons limited by these frequencies is highlighted with the green area in Supplementary Fig. 5). These phonons strongly affect the z component of the permittivity tensor, and thus they can couple with the electric field components in the z direction of the focused light, which explains the extra dip.
For a more detailed analysis of this coupling, we show in Supplementary Fig. 5d the calculated reflectivity spectra under focused illumination as a function of the total cavity thickness L cav . When the frequency of the bare-cavity mode gets close to the Reststrahlen band of the out-of-plane phonon, an anticrossing characterized by a relatively small splitting is observed, as indicated by the dots as a guide. Hence, the use of a focused beam makes possible to observe the coupling of the cavity modes with in-plane and out-of-plane phonons.
The coupling with the out-of-plane phonons has also been observed in the experimental spectra of the fully filled hBN cavities. Supplementary Figure 6a shows the measured reflectivity spectra for all six cavities. This figure corresponds to a zoom of Fig. 5c in the main text on the region near the Reststrahlen band associated with the out-of-plane phonons (green area). We observe a dip labelled by ω (1) cav that changes strongly with the thickness of the cavity and it is mostly associated to the bare-cavity mode. The other dip is close to ω LO, , and its frequency is also indicated by the yellow dots in Fig. 5a of the main text. In order to confirm the nature of the dips at frequency ≈ ω LO, , cav . On the other hand, once we consider focused illumination, we can observe the second feature near ω LO, that was identified in the experiments. We thus confirm that the extra peak only appears in the calculations when the polarization has a nonzero z component and thus the illumination couples also with the out-of-plane phonons. However, we notice that the size of the experimental dips is significantly larger than the simulated ones, and we attribute this discrepancy to experimental imperfections such as rugosities (which can scatter light at high angles) or non-perfect planarity of the fabricated cavities. Despite this difference, the good agreement on the spectral positions of the dips indicates that they are indeed a result of the coupling between the first cavity mode and the out-of-plane phonons.

Angle dependence of the reflectivity spectra
We show in Supplementary Fig. 7 the reflectivity spectra of hBN cavities under p-polarized light, as a function of the angle θ of the incident illumination with respect to the normal axis. Panel a corresponds to the cavity with thicknesses 510 nm/10 nm/370 nm of the MoS 2 /hBN/MoS 2 layers (see inset), which is analyzed in Fig. 2 in the main text. The results depend very weakly on the angle θ of the incident illumination because the corresponding angle of the refracted beam inside the cavity is much smaller than θ (due to the high permittivities of the materials filling the cavity). For reference, we show in Supplementary Fig. 7b the spectra of the bare cavity with the same thickness, which shows the same weak dependence on the angle.
Further, we show the spectra of a cavity of thickness L cav = 1665 nm fully filled with hBN ( Supplementary Fig. 7c) and with a material of permittivity ε hBN,∞ = 4.52 (Supplementary Fig.   7d). We observe again that the influence of the angle on the reflectivity is small. The main effect of increasing θ is the emergence in Supplementary Fig. 7c of new faint reflectivity dips close to the frequency of the out-of-plane hBN phonon ω LO, = 819 cm −1 . This additional dip indicates that for oblique incident angles it is possible to observe the coupling of the cavity with the out-of-plane phonon.
Last, we have checked that the reflectivity spectra under s-polarized light contains almost the same dips as in Supplementary Fig. 7. The main difference is that in this case there is no dip near ω LO, (i.e. it is not possible to couple with the out-of-plane phonon), because the electric field is polarized in the in-plane direction.
6 Analytical expression for the coupling strength 6

.1 Coupling strength of a TO phonon with a Fabry-Pérot cavity mode
In Supplementary Note 2.3, we describe the method that we follow in the main text to calculate the coupling strength g for hBN cavities, by extracting the polaritonic frequencies from the poles of the reflection and transmission coefficients and fitting these frequencies with a model of coupled harmonic oscillators. In this note, we present an alternative approach based on the microscopic interaction between the phononic material and the modes of a Fabry-Pérot cavity, which allows for deriving an approximate analytical expression for g, and can give additional physical insight about the coupling.
An analogous approach has been followed in Supplementary Ref. [18] to model the coupling strength between molecules and guided nanowire modes.
In this derivation, we consider a simplified picture where we model each unit cell of hBN as a harmonic oscillator characterized by a transition dipole moment d i . The dipole is induced by the vibrations of the atoms. From this perspective, the material is a collection of N cell dipoles interacting with the cavity mode. Due to the homogeneity of the material, the transition dipole moments of all unit cells are the same: d i = d. All dipoles are oriented parallel to the electric field, and in order to deduce the module d i = |d i |, we first consider the electric susceptibility of the material, which for a polar material such as hBN is ω TO and ω LO are the frequencies of the transverse optical and the longitudinal optical phonon, respectively, γ is the damping frequency and ε ∞ is the high-frequency permittivity. From this expression, we can directly obtain the polarizability α i of each dipole induced in a unit cell of volume V cell . We consider that the polarization density P is related to the electric field E as P = α i V cell E = ε 0 χE, and we therefore obtain where ε 0 is the vacuum permittivity. We can now relate the transition dipole moment d i and the classical polarizability by [17] and, from Supplementary Eqs. (17)- (19): This expression allows to characterize the coupling strength g i between the dipole associated to a particular unit cell and the cavity mode with quantized electric fieldÊ(r). This electric field corresponds to a bare cavity where the hBN is substituted by a material of frequency-independent, real permittivity ε ∞ (i.e., without considering the resonant polarizability of the unit cells). In this derivation, we consider normal incidence of the excitation light. For the quantization of the field, we further assume that the mirrors of the cavity are perfect, so that the electric and magnetic energy of the modes are equal and the fields do not penetrate into the mirrors. The electric field of the Fabry-Pérot mode of order j resonant at frequency ω (j) cav is quantized as [19] a andâ † are the annihilation and creation operators of the cavity mode, respectively. The fields only vary along the z direction normal to the flat interfaces of the system, as given by the field profile ρ(z) (normalized so that ρ(z) = 1 in the position where the field amplitude is maximum). The electric field is assumed to be polarized along the x axis, with unit vectorû x . V eff = S ε(z)|ρ(z)| 2 dz is the effective volume of the field, where S is the effective surface of the cavity and ε(z) indicates the spatial distribution of the permittivities of the system. The value of the electric field E(r) of a Fabry-Pérot mode is related to the coupling strength g i between that mode and the dipole excited in the unit cell at position r i , according to Taking into account that the dipole moments and the electric field are parallel, we write this coupling strength as We next consider the Hamiltonian describing the interaction of the cavity mode with the ensemble of dipoles associated with all the unit cells in the phononic material whereb i andb † i are the annihilation and creation operator of the harmonic oscillator that represents each vibrational dipole, respectively. The next step is to write the Hamiltonian of the system in terms of the annihilationb c and creationb † c operators of the bright collective excitation that represents the TO phonon,Ĥ = ω (j) where we have also added the annihilationĉ i and creationĉ † i operators of the N cell − 1 dark collective modes that are uncoupled to the cavity mode. Supplementary Eq. (25) can be derived from Supplementary Eq. (24) by using the Dicke transformation where g takes the value [20] where in the last step we have taken into account that the density of dipoles is dN cell dV = 1 V cell and the effective surface S of the mode is eliminated after the integration in the x − y plane. The integral of the denominator in Supplementary Eq. (28) is evaluated in the entire cavity, while the integral of the numerator is bounded to the region of the phononic material. Therefore, in this subsection we have obtained an expression that allows to calculate the coupling strength between the TO phonon of a polar material and a Fabry-Pérot mode of a cavity with an arbitrary spatial distribution of the permittivity ε(z) (and thus arbitrary field distribution ρ(z)).

Application to the hBN microcavities
In this subsection, we focus on the particular system that we analyze in the main text, a cavity of thickness L cav with MoS 2 as a spacer and containing a layer of hBN (details can be found in Supplementary Note 1). First, we need to calculate the field profile ρ(z) of the modes of the bare cavity, in order to evaluate the coupling strength g through Supplementary Eq. (28). Labelling the positions of the MoS 2 -hBN interfaces as L 1 and L 2 (indicated in Supplementary Fig. 1), the spatial distribution of the permittivity ε(z) of the bare cavity is described by the function In each interval with constant permittivity ε i , the field profile ρ(z) satisfies the Helmholtz equation where c is the speed of light in vacuum. Further, since the cavity is assumed to be bounded by perfect mirrors, the electric field vanishes at both ends: ρ(0) = ρ(L cav ) = 0. In order to verify the boundary conditions, ρ(z) needs to be continuous and differentiable in all interfaces, which leads to the solution A is the normalization constant chosen so that the maximum of the field profile is ρ(z) = 1, while ω (j) cav and φ (j) are the j th solution of the system of equations We can obtain g by inserting Supplementary Eq. (31) into Supplementary Eq. (28), with the integral extending over the phononic material (between L 1 and L 2 ), i.e. we evaluate In this calculation, we always assume that the cavity mode of order j is resonant with the TO phonon, cav , so that the cavity thickness L cav is different for each order and the permittivity of MoS 2 is evaluated at that frequency (ε MoS 2 ≡ ε MoS 2 (ω TO ) in this subsection). Adjusting the cavity parameters (specially the thickness of the hBN layer) allows for varying g from 0 to a maximum value g max corresponding to the fully filled cavity that can be evaluated to be g max = Interestingly, we observe that g max only depends on the frequencies of the transverse optical phonon and the longitudinal optical phonon. For hBN, g max = 428 cm −1 .
Further, we can also obtain a simple analytical expression of g for thin layers of hBN placed in the middle of the cavity and interacting with the first cavity mode. In order to obtain the field profile ρ(z) in this regime, we consider that the thin layer of hBN disturbs very weakly the electric field of the cavity fully filled with MoS 2 and with the same total thickness L cav . Under this assumption, Moreover, since the thin layer is located in the position of the maximum amplitude of the electric field, this amplitude varies slowly inside the hBN layer, and we can assume that the field is constant between L 1 and L 2 (ρ(z) ≈ 1). By evaluating the integrals in Supplementary Eq. (33) for L 1 = Lcav 2 − L hBN 2 and L 2 = Lcav 2 + L hBN 2 , we obtain This result is consistent with the Dicke model [21,22] that has been often applied to quantify the coupling strength between a cavity mode and an ensemble of N identical molecules [23][24][25]. Since in our case the thickness of the hBN layer L hBN is proportional to the amount of dipoles interacting with the cavity mode, we obtain the same dependency of g on the amount of matter excitations g ∝ √ N as predicted by Dicke.
We assess in Supplementary Fig. 8 the validity of our analytical model by comparing the result of Supplementary Eq. (33) with the coupling strengths calculated by applying the coupled harmonic oscillator model to the polaritonic frequencies calculated with the transfer-matrix method. For the calculation of the polaritonic frequencies, we always choose the thickness L cav so that the frequency of the cavity mode is resonant with the TO phonon frequency. We show the evolution of the coupling strength g as a function of the filling factor L hBN Lcav for the first three cavity modes. Generally, g increases with the filling factor, particularly strongly for thin layers and odd modes, until for all modes it saturates at the same value g max for fully filled cavities. The third mode also shows a nonintuitive tendency of the coupling strength which decreases with an increase of L hBN Lcav (in the interval 0.5 L hBN Lcav 0.7). The reason is that changing L hBN also modifies L cav (to keep the cavity resonant with the TO phonon) and the field distribution inside the cavity, which leads to a complex tendency. inside the cavity. The total thickness of the cavity is chosen so that the first cavity mode is resonant with the TO phonon frequency, and MoS 2 is chosen as the spacer material (see sketches). Displacing the 10 nm hBN layer from a position near a mirror (left of the x axis in the figure) to the center of the cavity (right of the axis) induces a transition from the weak to the strong coupling regime. In a similar way, the coupling strength for the cavity containing a 100 nm hBN layer is tuned along a range of values corresponding to the strong coupling regime. Last, for L hBN = 200 nm, the system is tuned from the strong to the ultrastrong coupling regime. Therefore, the results in Supplementary   Fig. 9a show that displacing a hBN layer of fixed thickness inside the cavity also makes possible to control the coupling strength [26].
The influence of the position of the hBN layer on the coupling strength is further emphasized in Supplementary Fig. 9b. We plot the coupling strength g as a function of the filling factor L hBN Lcav for different situations. We first show that, when the hBN layer is placed in the center of the cavity, the coupling strength with the second cavity mode (red line) is generally much smaller than with the first mode (blue line), because in this position the intensity of the cavity mode vanishes for the former and is maximum for the latter. Further, we consider a situation that maximizes the coupling strength with the second cavity mode (red dots). To obtain these results, we calculate the intensity