Floquet group theory and its application to selection rules in harmonic generation

Symmetry is one of the most generic and useful concepts in science, often leading to conservation laws and selection rules. Here we formulate a general group theory for dynamical symmetries (DSs) in time-periodic Floquet systems, and derive their correspondence to observable selection rules. We apply the theory to harmonic generation, deriving closed-form tables linking DSs of the driving laser and medium (gas, liquid, or solid) in (2+1)D and (3+1)D geometries to the allowed and forbidden harmonic orders and their polarizations. We identify symmetries, including time-reversal-based, reflection-based, and elliptical-based DSs, which lead to selection rules that are not explained by currently known conservation laws. We expect the theory to be useful for ultrafast high harmonic symmetry-breaking spectroscopy, as well as in various other systems such as Floquet topological insulators.


Supplementary note 1: DS constraints as eigenvalue problems
We show here how each DS constraint on a physical observable (derived in Eq. (22) in the main text) can be written as a set of eigenvalue problems in the Fourier domain. In order to understand how this takes place, we first demonstrate it for a general observable ( ) in a (2+1)D system with ̂2 =̂2 •̂2 DS. Since ( ) isperiodic ( -periodicity is a symmetry of all Floquet systems), it can be expanded to a Fourier series: (t) = ∑ n exp ( 2 ) (1) where the coefficients n are complex numbers. Thus, the Fourier transform of ⃗( ) is discrete in the frequency domain, with Dirac delta functions around integer multiples of the frequency 0 = 2 / . Furthermore, if ( ) is a real quantity, then, the coefficients real parts are even, and their imaginary parts are odd: n = n + n ; n = −n , n = − −n (2) where n and n are real. Applying the symmetry constraint on supplementary Eq. (1), that is, requiring that the observable is invariant to ̂2 , yields: ∑̂2 • n exp ( 2 ) = ∑ n exp ( 2 ) The coefficients of identical exponents in supplementary Eq. (3) must be equal. This reduces supplementary Eq. (3) to separate eigenvalue problems for each frequency component : ̂2 • n = n (4) The rotational operation ̂2 has two identical eigenvalues of value -1, with any eigenvector. The constraints in supplementary Eq. (4) are upheld if the temporal part of the symmetry (that translates by half a period) leads to an eigenvalue that corresponds to those of the matrix ̂2 , which only happens for odd values of . That is, the equation is solved by any eigenvector for odd values of , and has only a trivial solution for even values of . All DSs can be treated in a similar manner, while only the eigenvalues and eigenvectors that uphold the equations change for different DSs. These constraints restrict the spectrum of the physical observable, and therefore its temporal evolution as well. A solution is always guaranteed: either a non-trivial one, or a trivial one ( n = 0). For HG, a trivial solution results in "forbidden harmonics", that is, selection rules. More generally, selection rules can be expressed as any relation between the components of n , determined by the eigenvectors themselves. A similar derivation for a ̂3 symmetry which is present in counter-rotating bi-circular -2 laser fields, was presented in ref. 1.

Supplementary note 2: Selection rules in groups with multiple members
In this note we prove that the generators of a symmetry group are sufficient to describe the selection rules and other requirements in Floquet systems. This is intuitive, since the generators provide all the information about the group (other members in the group result in restrictions that are upheld by closure and associativity in the group). First, we consider a group with a single generator, ̂. The generator can be broken to spatial and temporal components: ̂=̂•̂ (5) The resulting eigenvalue problems that represent the symmetry constraints due to ̂ can be written as: The temporal part determines the eigenvalues, and the spatial part determines the operator to be diagonalized. A non-trivial solution exists only if the spatial and the temporal parts exactly match. Assuming a solution n that solves supplementary Eq. (6), symmetry operations that are powers of ̂ lead to eigenvalue problems with the same solutions. The 'th power leads to: where we used the commutativity of the spatial and temporal operators. Due to associativity in the group, the equation above is solved by the same n , with the eigenvalues of ̂ to the power of . Thus, symmetry constraints resulting from powers of ̂ coincide with those imposed by ̂.
Next, we consider two generators, ̂ and ̂ and their resulting symmetry constraints. Again, we assume a solution n that simultaneously solves the eigenvalue problems of ̂ and ̂ combined. In this case, the same solution also solves the eigenvalue problem defined by the symmetry operation that is ̂•̂: where we again used the commutativity of the spatial and temporal operators. Supplementary Eq. (8) can be thought of as a series multiplication, or as a single eigenvalue problem whose eigenvalues are multiplications of those of ̂ and ̂. Clearly, if all symmetry constraints due to powers of ̂ and ̂ and their products are upheld, then the same solution is valid for all symmetry constraints in the group, since the group is comprised of powers and products of its generators. At this point the theorem is trivially extended to any number of generators, and results in the following: given a dynamical group , one must only examine the selection rules arising from the generating symmetries.

Supplementary note 3: HG selection rules
We develop the HG selection rules for underlying DSs using the eigenvalue problem approach described in supplementary note 1 for the time-dependent polarization (Eq. (24) in the main text), which is even under timereversal (odd observables, e.g. currents, have slightly different selection rules). We begin with the (2+1)D case that describes DSs in collinear atomic/molecular HG, and then move on to the (3+1)D case that describes noncollinear and solid HG. HG selection rules in (2+1)D Time reversal symmetry -̂ Τ symmetry results in the following constraints: n = −n (9) A general solution demands the and components of n are in phase, meaning any driver invariant under time-reversal results in linearly-polarized only harmonics. Reflection-time reversal symmetry -̂ ̂ symmetry results in the following constraints: The selection rules dictate that the and components have a relative phase of π/2. This means any elliptical emission is allowed, as long as the polarization ellipse has a major/minor axis parallel to the reflection axis. Reflection-time translation symmetry -̂ ̂ symmetry results in the following constraints: ̂• n = n (11) Here we must separate to even and odd cases. The general solution dictates selection rules where even harmonics are linearly-polarized along the symmetry axis, and odd harmonics are linearly-polarized orthogonally to the symmetry axis. This selection rule was demonstrated experimentally for a cross-linear ω-2ω pumps driving HHG in refs. 2,3. Reflection-time translation followed by time reversal symmetry -̂ ̂ symmetry results in the following constraints: The resulting selection rules dictate that the and components have a relative phase of π/2. That is, any elliptical emission is allowed, as long as the polarization ellipse has a major/minor axis parallel to the reflection axis. The identical selection rules to the operator ̂ are not coincidental, and indicate that it does not matter which point in time upholds ̂ symmetry (the additional time-translation operator shifts the ̂ symmetry by /2).
Rotation by 180 0 -time translation symmetry -̂ ̂2 symmetry was previously found responsible for odd-only harmonic emission from one-color laser fields 4 , as we also show here. ̂2 symmetry results in the following constraints: ̂2 • n = n (13) The spatial rotation operator by 180 o degrees has two identical eigenvalues of value -1. A non-trivial solution then only exists for odd values of . This results in odd-only harmonics, compatible with previous derivations, with no particular limitation on the polarization. For instance, slightly elliptical fields uphold this symmetry and produce elliptically polarized odd harmonics 5,6 . Rotation by 180 0 -time reversal symmetry -̂ ̂ symmetry results in the following constraints: ̂2 • n = −n (14) The resulting selection rules dictate that the and components of n are in phase. Thus, any system upholding ̂ symmetry generates only linearly-polarized harmonics.
Rotation by 180 0 -time translation followed by time reversal symmetry -̂ ̂ symmetry results in the following constraints: ̂2 • n = −n (15) For which the resulting selection rules are identical to those imposed by ̂ symmetry, with similar arguments. Rotation-time translation symmetry -̂ ̂ symmetry results in the following constraints: The rotational operator has two eigenvalues: ± 2 / , with the eigenvectors: ( 1 ± ). The selection rules dictate only circularly polarized = ± 1 harmonics are emitted (for integer ), with alternating helicities. This derivation appeared in ref. 7.

Generalized rotation-time translation symmetry -̂,
This symmetry is exhibited in HHG driven by generalized two-color bi-circular EM fields 8,9 . ̂, symmetry results in the following constraints: The resulting selection rules dictate only circularly polarized = ± harmonics are emitted (for integer ), with alternating helicities, consistent with previous derivations from conservation laws 8,9 . Generalized elliptical rotation-time translation symmetry -̂, ̂, symmetry results in the following constraints: The generalized rotation-scaling operator has the same eigenvalues as the standard rotational operator, but different eigenvectors that depend on the ellipticity: ( 1 ± ). The resulting selection rules dictate that the emitted harmonics are elliptically polarized along the elliptical axis with alternating helicity, with the same underlying ellipticity , and only = ± harmonics are emitted (for integer ). HG Selection rules in (3+1)D The (3+1)D formulation is identical to the (2+1)D case, except that rotations must now revolve around a specific axis in 3D space (which we always assume is the z-axis), and reflections are no longer reflection axes, but reflection planes (which we always assume to be the xy plane). Additionally, n is now a 3D object. Time reversal symmetry -̂ Τ symmetry results in the following constraints: A general solution demands the , , and components of n are in phase, meaning any driver invariant under time-reversal results in linearly-polarized only harmonics, same as the (2+1)D case. Reflection-time reversal symmetry -̂ ̂ symmetry results in the following constraints: The selection rules dictate that the and components are in phase (linear polarization within the plane), and are phase shifted by π/2 for the -axis induced polarization. This means the 3D polarization ellipsoid has a 2D projection for each harmonic with a major/minor axis along the -axis. Reflection-time translation symmetry -̂ ̂ symmetry results in the following constraints: Here we must separate to even and odd cases. The general solution dictates selection rules where even harmonics are in any elliptical polarization within the symmetry plane, and odd harmonics must be linearlypolarized orthogonally to the symmetry plane. Reflection-time translation followed by time reversal symmetry -̂ ̂ symmetry results in the following constraints: The resulting selection rules are identical to those of the operator ̂. Rotation by 180 0 -time translation symmetry -̂ ̂2 symmetry results in the following constraints: ̂2 • n = n (23) The selection rules lead to odd-only harmonics that can be polarized within the plane orthogonal to the rotation ( plane), and only even harmonics polarized parallel to the symmetry axis ( -axis). Rotation by 180 0 -time reversal symmetry -̂ ̂ symmetry results in the following constraints: ̂2 • n = −n (24) The resulting selection rules dictate that the and polarization components are in phase (linear polarization within the plane), and phase shifted by /2 from the symmetry axis. Thus, the 3D polarization ellipsoid has a 2D projection for each harmonic with a major/minor axis along the -axis. Rotation by 180 0 -time translation followed by time reversal symmetry -̂ ̂ symmetry results in the following constraints: (25) For which the resulting selection rules are identical to those imposed by ̂ symmetry. Inversion-time translation symmetry -̂ ̂ symmetry results in the following constraints: (26) The selection rules lead to odd-only harmonics with any generalized polarization. We note this symmetry is well-known from perturbative nonlinear optics, and often attributed to centro-symmetric media which cancel out even order nonlinear coefficients 10,11 . In our case, ̂ symmetry is the generalization of this selection rule, which also depends on the DS of the incident pump field. Inversion-time reversal symmetry -̂ ̂ symmetry results in the following constraints: The resulting selection rules dictate that the , and axes induced polarization components are in phase (any linearly polarized high harmonic emission is allowed). Inversion-time translation followed by time reversal symmetry -̂ ̂ symmetry results in the following constraints: (28) For which the resulting selection rules are identical to those imposed by ̂ symmetry. Generalized elliptical rotation-time translation symmetry -̂, ̂, symmetry results in the following constraints: The resulting selection rules dictate that the emission of = ± (for integer ) harmonics is permitted with elliptical polarization and alternating helicity along the elliptical axis and within the plane orthogonal to the rotation axis, with the same underlying ellipticity . The emission of = harmonics is also permitted, but only with linear polarization parallel to the symmetry axis. This operation reduces to circular symmetry (̂, ) for = 1. Elliptical improper rotation-time translation symmetryeven orders -̂, ̂2 , symmetry results in the following constraints: The resulting selection rules dictate that the emission of = 2 ± (for integer ) harmonics is permitted with elliptical polarization and alternating helicity along the elliptical axis and within the plane orthogonal to the rotation axis, with the same underlying ellipticity . The emission of = (2 + 1) harmonics is also permitted, but only with linear polarization parallel to the symmetry axis. This operation reduces to circular improper rotational symmetry (̂2 , ) for = 1. Elliptical improper rotation-time translation symmetryodd orders -̂+ , ̂2 +1, symmetry results in the following constraints: The resulting selection rules dictate that the emission of = 2 (2 + 1) ± (for integer ) harmonics is permitted with elliptical polarization and alternating helicity along the elliptical axis, and within the plane orthogonal to the rotation axis, with the same underlying ellipticity . The emission of = (2 + 1)(2 + 1) harmonics is also permitted, but only with linear polarization parallel to the symmetry's rotation axis. This operation reduces to circular improper rotational symmetry (̂2 +1, ) for = 1.

Connection between DS group theory and perturbative selection rules
In this subsection we directly outline the connections between our general DS-based selection rules derivation for HG, and the standard perturbative derivation of selection rules in non-linear optics (NLO). We also discuss the connection between the medium's static symmetries and the selection rules presented in the main text.
First, we comment on the role of the medium's characteristic point-group symmetry (describing the static symmetries of the field-free system, either gas, liquid, or solid) in the DS-based theory. In practice, the combination of the crystal/molecule symmetry and pump symmetry is done by taking the following steps: (i) Identify the symmetry of the medium using standard tools (molecular group theory for molecules, and crystallographic space group theory for solids). (ii) Identify the DSs of the pump. (iii) Identify the DSs of the joint system (pump and medium), that is, the static symmetries of the medium that are also elements in the DSs of the pump. (iv) Find the selection rules corresponding to the DSs of the joint system in the tables presented in the main text.
This approach is general, and can be applied to any target with one of the 32 known point groups systematically according to the symmetry element. As an example, consider an α-Al 2 O 3 single crystal interacting with an -2 bicircular pump 12 propagating along the c-axis of the crystal. In this case, we identify α-Al 2 O 3 belonging to the R3c space group that exhibits 3-fold symmetry (̂3) along the c-axis 13 . Similarly, we identify the bi-circular pump's dynamical symmetry group that exhibits ̂3 DS, which have an identical ̂3 spatial part. Hence, ̂3 is a joint DS for the light-matter system, leading to HG selection rules according to Table 1 in the main text. On the other hand, a single-crystal of rocksalt MgO belongs to the space group Fm3m that does not exhibit ̂3 symmetry; hence, if MgO is pumped by the same bi-circular beam, ̂3 is not a joint DS of the light-matter system, and does not lead to selection rules.
Second, we comment on the general theoretical connection between the perturbative and DS-based approaches: standardly in NLO, the symmetries of the medium (either gas, liquid, or solid, which are known from the specific point-group of the medium) are used in order to reduce the non-linear optical coefficients tensor to its minimal representation (to the smallest number of nonzero independent parameters). Using this approach, selection rules for specific perturbative NLO processes can be derived based on the elements of the reduced nonlinear optical tensor 10,11 . Our DS group theory based derivation generalizes this approach, giving rise to a unified description of selection rules for both perturbative NLO processes, and the non-perturbative processes (e.g., HHG).
Third, we discuss generalization of several selection rules by our approach: (i) Generation of even-order harmonics is forbidden in centro-symmetric media 10,11 . This selection rule is re-derived for Hamiltonians exhibiting the DS ̂. More generally, we find that the same selection rule is derived for Hamiltonians exhibiting the DS ̂2 , such that even if the medium is not centrosymmetric and only exhibits ̂2 symmetry along the pump's propagation axis, the same selection rule is observed. (ii) Rotationally-invariant crystals pumped by a monochromatic circularly polarized field generate only circularly polarized harmonics 14 . This selection rule is re-derived for Hamiltonians exhibiting the DS ̂ (even when the pump field is not monochromatic). We also generalize this selection rule to the selection rules derived for Hamiltonians exhibiting the DSs ̂, ,̂, . (iii) No HG in a spherically-symmetric media pumped by a monochromatic circularly polarized field 14 .
This selection rule is re-derived for Hamiltonians exhibiting the DS ̂ for →∞. This means that this selection rule is also observed from media exhibiting a continuous rotational symmetry along the pump's propagation axis (cylindrical symmetry).

Supplementary note 4: Examples of dynamical groups in (2+1)D
Examples of dynamical groups in (2+1)D are given in supplementary

Supplementary note 5: Similarities to other groups
As mentioned in the main text, dynamical Floquet groups are similar to those of line-groups 15 for the (2+1)D case. Considering line-groups that extend an infinite and periodic -axis, this axis is similar to a periodic and infinite time-axis. Translations along the -axis in line-groups are equivalent to time-translations (̂). All operations within the xy plane are identical to the spatial symmetry elements in (2+1)D dynamical groups. The main difference is in the time-reversal operation (Τ), which geometrically is similar to a spatial reflection plane within the plane in line-groups. However, time-reversal is an anti-unitary operation unlike reflections. We further note that certain features of space-time groups 16 , knot-groups 17 , and space-groups 18 , are isomorphic to Floquet groups. Due to the resemblence, we do not derive here the full list of character tables and multipcation tables of Floquet groups, but refer the reader to refs. [15][16][17][18] , and references therein.

Numerical model
We numerically solve TDSE (2D/3D) in the length gauge, within the single active electron approximation, and the dipole approximation. In atomic units the TDSE is given by: where | ( )⟩ is the time dependent wave function of the single electron and ( ) is the atomic potential. We use a spherical atomic model, where the electron is initially in the 1s orbital ground state found by complex time propagation. The Coulomb interaction is softened at the origin with two free parameters set to describe the ionization potential ( ) of the Ne atom ( = 0.793 a.u.), and He + ion ( = 1.999 a.u.). We set the atomic potential according to ref. 19: where = 1, 2, and = 0.1195, 0.1583 a.u., for Ne and He + in the 2D case, respectively, and Z= 1.285, a=0.0001 for Ne in the 3D case. To avoid reflections from the boundaries, we use a complex absorbing potential set to:  Fig. 1 for several exemplary cases. Each harmonic spectra shows the and components of the emitted harmonics with the numerically calculated ellipticity at each peak. The different drivers are presented in the inset of each harmonic spectra as a parametric Lissajous curve. The components of the field are given up to an overall amplitude coefficient. The analytically derived selection rules for drivers with ̂ symmetry are such that all harmonics are linearly-polarized, where even and odd harmonics are oriented along, and perpendicular, to the symmetry axis, respectively. As seen in supplementary Fig. 1, the harmonic selection rules are generally upheld very well by all harmonic orders, yet deviations can appear near resonances, which are marked with dashed lines. For instance, as shown in supplementary Fig. 1(c), the ellipticity of the 10'th harmonic, which is near a resonance, is 0.09 instead of 0. As shown in supplementary Fig. 2, the deviation reduces when the length of the pulse or the length of the turnon period of the flat-top envelope are increased. Thus, we conclude that resonances increase the sensitivity of HHG to broken time translation symmetry (as exists in finite pulses)an important and potentially useful feature in HHG spectroscopy through DS breaking. Supplementary Fig. 1 HG spectra for driver fields with ̂ DS from a model Ne atom system. a-f HG spectra from various types of driver fields that exhibit ̂ DS. The driver field is indicated in inset up to an overall amplitude coefficient with its Lissajous curve. Atomic resonances are indicated by dashed lines, where the last line indicates the Ip. According to the selection rules derived from group theory, all harmonics should be linearly polarized, even harmonics along the symmetry axis, and odd harmonics orthogonal to it. The calculated ellipticity is presented at each peak, indicated in black, and intensity in the y-axis is given in log scale.  Fig. 1(c) is shown as a function of changes in the pulse envelope that perturb time-translation symmetry. a Ellipticity as a function of the total pulse duration for an envelope with a rise and drop sections 10 optical cycles long, and a flat top of varying duration. b Ellipticity as a function of varying turn-on duration in the flat-top envelope, with a flat-top part that is 10 optical cycles long. With increasing pulse duration (and turn-on duration), the ellipticity indeed approaches zero, as predicted by the selection rules for the symmetry in this system.

Time-reversal related symmetries (̂,̂,̂) Next, we examine (2+1)D DS operators that involve time-reversal symmetries (Τ ̂,̂= Τ •̂2,̂= Τ •̂2 •̂2).
According to group theory, these DSs should result in only linearly-polarized harmonics. However, these symmetries are slightly broken due to tunneling ionization (see main text). We numerically examine two systems with different ionization potentials: Ne and He + . In He + , the electron is tightly bound by the deep potential well, which results in a high , and a large number of below harmonics that do not break timereversal symmetry (since ionization does not play a role in their generation). Various examples are presented in supplementary Figs. 3 and 4. Each harmonic spectra shows the harmonic intensity and the numerically calculated ellipticity at each peak. The different drivers are presented in the inset of each harmonic spectra as a parametric Lissajous curve. The -components of the field are given up to an overall amplitude coefficient. As shown in supplementary Figs. 3 and 4, the harmonic selection rules are generally upheld by harmonic orders below the , where deviations are mostly found near resonances. Above , the predicted selection rules are sometimes upheld and sometimes not. For example, in the spectrum of supplementary Fig. 3(b), several high harmonics have a relatively large ellipticity (up to 0.5) due to the broken symmetry. In contrast, in the spectrum of supplementary Fig. 4(b), all harmonics are practically linearly-polarized. This broken symmetry thus depends on several parameters, such as the shape and intensity of the driver and of the atomic potential. Supplementary Fig. 3 HG spectrum for driver fields with DSs that involve time-reversal from a model Ne atom system. a a driving field with ̂ DS, b a driving field with ̂ DS, c and d different driving fields with ̂ DS. The DS is indicated in the inset of each spectrum with the driver field up to an overall amplitude coefficient, and its Lissajous curve. Atomic resonances are indicated by dashed lines, where the last line indicates the Ip. According to the selection rules derived from group theory all harmonics should be linearly-polarized. The calculated ellipticity is presented at each peak, indicated in black, and intensity in the y-axis is given in log scale.

th order elliptical symmetry (̂)
The kinetic energy and atomic/molecular potential terms are generally not invariant under the elliptical symmetry operation, ̂. Still, there are cases where a driver pump can be chosen to uphold an elliptical symmetry, and not break the spherical symmetry of other terms in the Hamiltonian. This can only happen for n=4, since every ellipse has four points that are equidistant from the origin (supplementary Fig. 5). For example, in ref. 23 a train of linearly-polarized pulses was proposed that upholds this symmetry relation, and it was indeed shown numerically to yield the predicted selection rules. In this case, within each pulse the electron's motion, and therefore the polarization, is confined to the polarization axis of the pulse. Thus, the polarization is invariant under the elliptical DS. Below, we present an intriguing case that also yields the predicted selection rules of the (2+1)D elliptical DS, even though the electron's motion is two dimensional and intricate. We explore HHG driven by a synthetic field defined by a piecewise function that upholds 4'th order elliptical symmetry. The driver field has the form: where 0 is the field amplitude, ( ) is a flat-top envelope with a three optical cycle ramp up and down and a 15 optical cycle long flat-top, is the ellipticity of the underlying symmetry (from 0 to 1), and 0 = 2 / is the optical frequency. Notably, producing this field requires an infinitely broad spectrum due to the non-smooth orientation shifts (the field is continuous but its derivative is not). Still, while such a field is not accessible experimentally, it can be approached by using multi-color drivers. This field provides 4 th order elliptical symmetry present in all terms in the Hamiltonian. This example is slightly more general than the case of linearly- calculated ellipticity, and the (+) and (-) projected helical components to indicate the helicity of the harmonics. The Lissajou parametric curve in each case is shown in inset, along with the ellipticity of the driver's underlying symmetry. In all cases the selection rules derived from group theory match with the calculated spectrum, except near atomic resonances. Supplumentary Fig. 5 Schematic illustration of 4 th order elliptical DS (̂) which still maintains circular DS for any ellipticity. In every ellipse four points equidistant from the origin exist, which belong both to a circle and an ellipse simultaneously. The red arrows represent scaling transformation for the y axis from the blue unit circle to the yellow ellipse. Supplementary Fig. 6 HG spectrum from a 4 th order elliptical DS piecewise electric driver field, as in supplementary Eq. (37). a-d, underlying ellipticity ε0=1,0.7,0.3,0.1, respectively. The ellipticity of the underlying symmetry is indicated in inset, as well as the field Lissajous curve. Calculated ellipticity for each harmonic peak is indicated in black, and according to the selection rules derived from group theory should be identical to the ellipticity of the driver's symmetry, with alternating helicities. Atomic resonances are indicated in dashed lines, where the last line indicates the Ip. Intensity in the y-axis is projected onto the left and right circularly polarized components, in red and blue, respectively, and the intensity is given in log scale.

(3+1)D Rotational symmetry in non-collinear HG (̂, )
Next, we examine (3+1)D rotational DS operators (̂, =̂•̂, ). These differ from the (2+1)D case by a specific choice of rotational symmetry axis, which we denote the -axis. According to group theory, this DS should result in counter rotating ± harmonics that are circularly polarized within the xy plane, and also harmonics that can only be emitted with linear polarization along the -axis (for integer ). Fascinatingly, emission of the harmonics means that there is photon mixing between photons with all possible polarization axes (since otherwise the even harmonics would be forbidden). We numerically verify this in a Ne like 3D system for a non-collinear geometry with ̂3 symmetry, where HG is driven by a counter rotating ω-2ω bicircular pump propagating along the z-axis, and an orthogonally propagating (along the x-axis) pulse of frequency 3ω polarized linearly along the z-axis, seen in supplementary Fig. 7. The spectrum shows the harmonic intensity projected onto the circular components within the xy plane with the numerically calculated ellipticity at each peak, as well as the remaining -axis emission. The driving field is presented in the inset as a parametric Lissajous curve, and schematic illustration of the non-collinear setup.
Supplementary Fig. 7 HG spectrum from a (3+1)D ̂ symmetric pump. The driving field is given in the inset up to an over amplitude coefficient. Red and blue stand for the circular projection within the xy plane, and calculated ellipticity for each harmonic peak within is indicated in black, and according to the selection rules should be exactly 1, with alternating helicities. The z-axis emission is presented in green, and is shown to be only at the predicted harmonic frequencies 3q (for integer q). Intensity is given in log scale.

(3+1)D Improper Rotational symmetry in noncollinear HG (̂, )
Next, we examine (3+1)D rotational DS operators (̂2 , =̂2 •̂2 , •̂ℎ). According to the presented theory, this DS results in counter rotating 2 ± harmonics that are circularly polarized within the xy plane, and also (2 + 1) harmonics that can only be emitted with linear polarization along the -axis (for integer ). The emission of the (2 ) harmonics is symmetry forbidden, which is very surprising because this gives the appearance that selection rules from each axis are upheld separately, as if no mixing occurs. As an example, we explore ̂4 ,1 symmetry in a non-collinear geometry, with a bi-circular -3 counter-rotating pump pulse propagating along the -axis, and an orthogonally propagating pump pulse (propagating along the x-axis) that is linearly polarized along the -axis at frequency 2 . In this case, harmonics 3,5,7,9,11… should be circularly polarized with alternating helicities within the plane. Harmonics 2,6,10,14,18… should be emitted only with polarization components along the -axis. The remaining 4,8,12,16… harmonics are symmetry forbidden. We numerically verify this in a Ne like system, seen in supplementary Fig. 8. The spectrum shows the harmonic intensity projected onto the circular components within the plane with the numerically calculated ellipticity at each peak, as well as the remaining -axis emission. The driving field is presented in the inset as a parametric Lissajou curve, and schematic illustration of the non-collinear setup.  Supplementary Fig. 8 HG spectrum from a (3+1)D ̂, symmetric pump. The driving field is given in the inset up to an over amplitude coefficient. Red and blue stand for the circular projection within the xy plane, and calculated ellipticity for each harmonic peak within is indicated in black, and according to the selection rules should be exactly 1, with alternating helicities. The z-axis emission is presented in green, and is shown to be only at the predicted harmonic frequencies 2(2q+1) for integer q. Intensity is given in log scale.

Supplementary note 7: HG selection rules from molecules and DS breaking spectroscopy
Here we present numerical examples for molecular HG selection rules, and for molecular symmetry and orientation spectroscopy from DS breaking. Specifically, we consider a planar 3-fold symmetric molecule that belongs to the 3ℎ point group (which supports 120 0 rotations, 180 0 rotations, and also reflections along the major molecular axes, within the BO approximation), interacting with a pump field that exhibits ̂=̂2̂ DS. We demonstrate that the selection rules derived from Floquet group theory depend both on the pump's DS, and on the symmetries of the molecule when it is oriented, and that the relative alignment with respect to the pump field determines whether a selection rule is observed (which is used for orientation spectroscopy). In case the medium is non-oriented, we show that the symmetry group of the orientation averaged ensemble and the pump's DS determine the observed selection rules.
According to the Floquet group theory, we should observe ̂ DS based selection rules in the HG spectrum from oriented molecules if the reflection axis of the pump is co-aligned with that of the molecules, leading to linearly polarized harmonics, where even/odd harmonics are polarized parallel/perpendicular to the symmetry axis. When the axes are not aligned, symmetry breaking will occur due to the molecular potential, and the selection rule will be broken. For the random non-oriented medium, theory predicts that the selection rule will be observed, because the orientation averaged ensemble is O(3) spherically symmetric (as discussed in the main text, and will only differ for chiral media), and supports the required reflection symmetry.
We consider the following 2D molecular potential that describes a poly-atomic molecule: where we set ( , ) = 0 (cos( + ) , sin ( + )), for 0 = 2.5 a.u., = {0,120 0 , 240 0 }, i=1,2,3, and we set as the relative angle of the molecule with respect to the -axis. This potential is represented on a realspace Cartesian × grid for = 120 a.u., with grid spacing Δ = Δ = 0.1171 a.u., where the ground state and first excited states are found by complex time propagation. The TDSE is solved as in supplementary note 6 for the case where the molecule is interacting with the following cross-linear -2 driving field that possesses the reflectional ̂ DS, but not 3-fold rotational DSs: ( ) = 0 ( )(sin( 0 )̂+ sin(2 0 )̂) where ( ) is a flat-top envelope with a 4 optical cycle ramp up and down and an 6 optical cycle long flat-top, 0 = 2 / is the optical frequency that corresponds to an 800nm wavelength, and 0 is chosen such that the maximal intensity of the pulse is 10 13 [W/cm 2 ]. The TDSE is solved for many relative orientations with respect to the laser field, such that we obtain the HG spectrum as a function of and the HG spectrum from the random non-oriented media (see supplementary Fig. 9). Supplementary Fig. 9(a) shows that when the molecular axis is aligned with the pump's reflection axis, the selection rules for ̂ symmetry are observed. On the other hand, these selection rules are broken if the axes are not aligned (supplementary Fig. 9(b)). Furthermore, for the non-oriented medium, the selection rules from the pump's symmetry are upheld regardless of the molecular medium (supplementary Fig. 9(c)). In fact, most HG experiments deal with such non-oriented media, meaning usually the observed selection rules are identical to those observed from atomic media 24 (this is why the majority of numerical examples above are presented for the atomic case). We also note that the molecule's symmetry properties that are not exhibited by the pump field are washed out in the spectrum, that is, there is no selection rule due to 3-fold symmetry even if the medium is aligned. This is because a symmetry based selection rule is only observed if the full Hamiltonian supports the DS, and in this case the pump is not 3-fold symmetric. Importantly, similar results are obtained when the TDSE is solved with other initial states (and also from degenerate states, as long as an equal population of degenerate states is chosen), and in a full 3D molecular potential. All of these results match the predictions by the Floquet group theory.
Next, we consider DS breaking spectroscopy for the molecule with ̂ symmetry -we can identify the molecular orientation in aligned media by examining when the selection rule from the pump is restored (orientation spectroscopy was demonstrated experimentally in SF 6 molecular media through a rotational DS in ref. 24). This is shown in supplementary Fig. 10(a) for the average intensity of the even/odd harmonic spectrum polarized along the / axes, respectively, and in supplementary Fig. 10(b) for the average ellipticity of the spectrum (harmonics 3-15). As seen, when the molecule is aligned with the pump, the selection rule is restored (hence can be measured), and we can also identify the presence of three reflection planes in the molecule. Conceptually, one can employ multiple pump fields of varying DSs and deduce the full molecular point group in this manner. Supplementary Fig. 9 Molecular HG spectrum from a ̂ symmetric pump and 3-fold molecular potential. The driving field's Lissajou plot is given in the inset together with the molecular orientation. The calculated ellipticity for each harmonic peak is indicated in black. Intensity is given in log scale. a Symmetric case, b symmetry broken case, c orientation averaged spectra (symmetric case).
Supplementary Fig. 10 Molecular HG DS breaking orientation spectroscopy from a ̂ symmetric pump and a 3-fold molecular potential. The driving field's Lissajou plot is given in the inset together with the molecule's orientation. a Intensity spectroscopythe selection rules from the pump field dictate that x-polarized harmonics are only odd, and y-polarized harmonics are only even, which is upheld when the molecule's reflection axis is aligned with the pump. b Ellipticity spectroscopythe selection rules from the pump field dictate that the spectrum is fully linearly polarized (0 ellipticity), which is upheld when the molecular axis is aligned with the pump.

Supplementray note 8: HG selection rules from solids and DS breaking orientation spectroscopy
Here we give some examples for new selection rules for HG from solids. In this particular case, the symmetries of the solid and the respective angle of incidence of the pump can determine the symmetries of the system. For example, a monochromatic laser field of any ellipticity upholds the DS ̂=̂2 •. If the solid has an inversion center (centro-symmetry) is invariant, and generation of even harmonics become forbidden. In fact, this is a well-known symmetry based selection rule in nonlinear optics 10,14 that is usually attributed to the symmetry of the solid, and our derivation presents its generalization. A second noteworthy example is that of the DS ̂=̂2 •̂, which is a symmetry of any monochromatic linearly polarized laser field. When driving a solid with a reflection symmetry plane that is aligned orthogonal to the pump's polarization axis is again invariant, leading to orthogonally linearly polarized even and odd harmonics. Importantly, if the solid is tilted out of plane with the incident wave, symmetry breaking occurs, which allows spectroscopically finding the solid's orientation. Lastly, DSs could be used to probe spin-orbit coupling or magnetic interactions in the solid, since these terms could break certain DSs.
As a numerical example, we show that the solid's orientation can break the reflection DS ̂ if it is not aligned to the pump's orientation, and use this to demonstrate orientation spectroscopy in solid HG. We calculate the harmonic spectrum from a single electron in a model 2D Mathieu type lattice potential 25 ( , ) = − 0 (1 + cos (  2 )) (1 + cos (  2 )) where 0 = 0.25, = 10 a.u. This is done similarly to the atomic systems, on a real-space Cartesian grid spanning 20 × 20 lattice sites of size × for = 370 a.u., with grid spacing Δ = Δ = 0.2645 a.u. The ground state is found by complex time propagation that corresponds to a real-space delocalized wave function that populates the first valence band, as is developed in ref. 26. The TDSE is then solved as in supplementary note 6 with the following cross-linear -2 driving field that possesses only ̂ DS: ( ) = 0 ( )(sin( 0 + )̂+ sin(2 0 )̂) (41) where = /3, ( ) is a flat-top envelope with a 5 optical cycle ramp up and down and an 8 optical cycle long flat-top, 0 = 2 / is the optical frequency that corresponds to an 800nm wavelength, and 0 is chosen such that the maximal intensity of the pulse is 5 × 10 13 [W/cm 2 ]. The harmonic spectrum is calculated and seen in supplementary Fig. 11 projected onto and polarization components for two cases: (a) the solid is aligned with the pump and the two share a -axis reflection DS, and (b) the solid is rotated by 25 0 with respect to the DS axis (this Mathieu potential supports three reflection axes). From supplementary Fig. 11, clearly when the DS is upheld the spectrum contains orthogonally polarized even and odd harmonics only, but for the broken symmetry different polarizations are emitted for different harmonics. We use this for orientation spectroscopy in supplementary Fig. 11(c), which shows the calculated average ellipticity of harmonics 10-18 in the spectrum as a function of the solid's angle with respect to the pump field. The average ellipticity has strong minima at < + + = 0 0 , 45 0 , 90 0 , which are the exact angles of reflection planes in the solid. Hence, an experimental ellipticity measurement may deduce the solid's orientation through DS breaking spectroscopy, and also the existence or lack of reflection planes in the solid. Notably we also see a peak in the average ellipticity at ≈ 25 0 , which is the orientation farthest away from any reflection axis, meaning that the maximal symmetry breaking may also be a useful observable.