Symmetry Breaking and Spatiotemporal Pattern Formation in Photonic Time Crystals
Egor I. Kiselev Yiming Pan



Physics Department, Technion, 320003 Haifa, Israel
The Helen Diller Quantum Center, Technion, Haifa 3200003, Israel 
kiselev.physics@gmail.com


School of Physical Science and Technology and Center for Transformative
Science, ShanghaiTech University, Shanghai 200031, China

In this work, we explore the dynamics of time varying photonic media
with an optical Kerr nonlinearity and an associated phase transition.
The interplay between a periodically modulated permittivity and the
nonlinearity induces a continuous transition of electromagnetic waves
to a state with broken spatial and time translation symmetries. This
transition gives rise to a lattice-like wave pattern, in many ways
similar to a spatial crystallization in solids. Symmetry breaking
triggers the emergence of soft, Goldstone-like modes, which propagate
as deformations of the lattice structure, as well as massive Higgs-like
modes – spatially uniform oscillations of the field amplitude. We
extend the analysis of the non-equlibrium symmetry breaking to 2+1
dimensional time varying media and discuss pattern formation as well
as the connection to discrete time crystals.

Symmetry Breaking and Spatiotemporal Pattern Formation in Photonic
Time Crystals
    Yiming Pan
    April 29, 2024
================================================================================





§ INTRODUCTION


The propagation of electromagnetic fields through media with time
varying material parameters has been considered early on <cit.>.
More recently, time varying media have attracted attention in photonics
<cit.>, metamaterials <cit.>,
and material science <cit.>, particularly
in the field of epsilon near zero materials <cit.>.
They are viewed as promising candidates for realizing such effects
as momentum-gapped (k-gapped) states, nonreciprocity <cit.>,
time switching <cit.>, time-varying
mirrors <cit.>, lasing and amplification
<cit.>. From a fundamental
perspective, the non-trivial statistical properties of randomly driven
time-varying media <cit.>,
their topological <cit.> and
radiative <cit.>
properties have been explored. Beyond optics, hydrodynamic <cit.>
and acoustic <cit.> time varying systems
have been studied. 

In this paper, we study nonlinear photonic time crystals (nonlinear
PTCs) – dielectric media with a periodically modulated permittivity
<cit.>. It is crucial
to distinguish PTCs from “time crystals” in the sense of Wilczek
<cit.>.
While it is important to maintain this distinction, we demonstrate
that, in the presence of a Kerr term, a PTC evolves to a subharmonic
steady state oscillating at half the driving frequency, and thus breaks
the discrete time translation symmetry of the drive on long timescales.
The steady state is very much akin to Faraday waves <cit.>,
which are in many ways similar to classical discrete time crystals
<cit.>.

One of the most striking features of PTCs is the so called momentum
gap (sometimes referred to as k-gap or q-gap, see Fig. <ref>a)
– an interval of wavenumbers in which the real part of the dispersion
ω(q) becomes flat and the imaginary part changes
sign, leading to exponentially growing modes inside the gap <cit.>.
In nonlinear PTCs and time varying media, the k-gap can lead to intriguing
effects such as the formation of superluminal solitons <cit.>.
We show that in PTCs featuring a generic Kerr nonlinearity, this growth
heralds an instability, leading to the emergence of non-equilibrium,
symmetry breaking steady-states: When the PTC's permittivity is modulated
at a frequency 2Ω, the field adopts a standing wave lattice
pattern whose wavenumber q^* and lattice constant λ=2π/q^*
are determined by the condition for parametric resonance q^*=Ω/c
(Fig. <ref>).

Furthermore, we show that in the symmetry breaking state, slow, long-wavelength
electromagnetic fields propagate through the PTC as distortions of
the lattice, locally contracting or expanding the lattice constant
λ (see Fig. <ref>). This behavior
is reminiscent of phonons in a crystal lattice, and the corresponding
distortions can be thought of as the Goldstone modes of the symmetry
broken state. Additionally, the symmetry broken state hosts spatially
uniform oscillations of the field amplitude at a characteristic frequency
that depends on the driving strength (see Fig. <ref>b).
These oscillations are gapped, massive modes that resemble, e.g.,
the Higgs modes of a superconductor.



§ RESULTS




 §.§ Symmetry breaking transition.




In this paragraph, we investigate the symmetry breaking transition
of nonlinear PTCs. We show that the transition is continuous when
the modulation strength of the dielectric constant exceeds a threshold
value determined by the damping.

The propagation of electromagnetic fields in a nonmagnetic, isotropic,
one-dimensional PTC with a generic, repulsive Kerr-nonlinearity is
described by the equation <cit.>

    1/c^2∂^2D/∂ t^2+γ/c^2∂ D/∂ t-(1+hcos2Ω t)∂^2D/∂ x^2=β|D|^2∂^2D/∂ x^2,

where D(x,t) is the displacement field, γ is
the decay rate of the electromagnetic field, c is the speed of
light and h is a dimensionless parameter controlling the modulation
strength of the dielectric constant. In what follows we assume h>0.
For brevity, we omit vector signs. For a small h≪1, an approximate
solution to Eq. (<ref>) can be found with the slowly varying
envelope approximation. Anticipating that the time dependent term
will lead to a parametric response at half the modulation frequency
Ω <cit.>, we choose the ansatz

    D=a(t)cos(Ω t)cos(q^*x)+b(t)sin(Ω t)cos(q^*x).

Here q^* is a critical wavenumber determined by the condition

    Ω=c|q^*|.


To first order in h, and projecting the nonlinear term of Eq. (<ref>)
onto the ansatz of Eq. (<ref>), i.e. disregarding
harmonics oscillating at frequencies ±3Ω and faster, we
obtain 

    ȧ    =-γ/2a-h̃Ω b+β̃Ω b(a^2+b^2)
    ḃ    =-γ/2b-h̃Ω a-β̃Ω a(a^2+b^2),

with h̃=h/4 and β̃=9β/32.

We deduce from the amplitude equations (<ref>) that
for h<h_c, where 

    h_c=2γ/Ω,

the PTC is stable around the trivial fixed point D(x,t)=0.
For h>h_c this fixed point becomes unstable, leading an initially
exponential growth of the amplitudes a, b. The nonlinear terms
then become increasingly important, and force the amplitudes to saturate.
Eventually, the PTC reaches a new fixed point, which corresponds to
the symmetry breaking phase that we want to study. 

It is useful to consider two distinct regimes. In the first regime,
which we will call this the near-critical regime, h is
close to the critical threshold value h_c. In this case, both,
the transition speed and the amplitudes a, b at the fixed point
are determined by the small parameter

    ϵ=h-h_c/h_c.

In the second regime, h≫ h_c holds. Here, the damping can be
disregarded to a good approximation. Let us term this the weakly-damped
regime. With γ set to zero, the equations (<ref>)
can be derived from the Hamiltonian 

    H(a,b)=h̃Ω/2(a^2-b^2)+β̃Ω/4(a^4+b^4)+β̃Ω/2a^2b^2.

Formally, we can think of a as a momentum variable and b as
a coordinate. The phase profile corresponding to H(a,b)
and the trajectories of a and b during the symmetry breaking
phase transition are shown in Fig. <ref>c.
The minima of the Hamiltonian (<ref>) are the
stable fixed points of Eqs. (<ref>), and are located
at

    a_0=0, b_0=±√(h̃/β̃).

For a small damping γ≪ hΩ, the solutions (<ref>)
are slightly modified. To first order in γ/hΩ, we find
a_0≈∓γ/Ω√(8β̃h̃)
and b_0≈±√(h̃/β̃)∓γ^2/16Ω^2√(β̃h̃).
For the most part (except when studying the transition to the symmetry
breaking state), we will not be interested in the damping, and only
include it for the sake of numerical stability. For weak damping,
the system will evolve to a state described by the above solution,
and the D-field will be given by

    D_0(x,t)=a_0cos(Ω t)cos(q^*x)+b_0sin(Ω t)cos(q^*x).

This steady state of the nonlinear PTC breaks the continuous spatial
translation, and the discrete time translation symmetries of Eq. (<ref>).
The latter follows from the fact that the solution D_0(x,t)
of Eq. (<ref>) is, unlike Eq. (<ref>) not invariant
under the time translation t→ t+π/Ω.

Let us now turn to the weakly damped regime where the modulation strength
h is close to the instability threshold h_c. This regime is
useful to study the transition dynamics from the trivial to the symmetry
breaking state. In particular, it can be shown that the transition
into the symmetry breaking steady state is a continuous phase transition.
To this purpose, we rewrite the amplitude equations for a(t),
b(t) in therm so the new variables a̅(t)=a(t)+b(t),
b̅(t)=a(t)-b(t). A perturbative
expansion of Eqs. (<ref>) in the small parameter ϵ
of Eq. (<ref>) (see supplement) then predicts that
the transitional dynamics can be described by a gradient descent equation
for b̅(t):

    ḃ̅̇=-∂ V(b̅)/∂b̅,

where V(b̅) is an effective double-well potential:

    V(b̅)=-ϵh̃_cΩ/2b̅^2+β̃^2Ω/48h̃_̃c̃b̅^6.

We plot the effective potential for different ϵ in Fig.
<ref>b and conclude that the transition into the
symmetry breaking state is continuous, i.e. the system goes through
a soft bifurcation when h reaches the critical driving strength
h_c. In terms of the original a and b, the potential minima
(and approximate fixed points of Eqs. (<ref>)) are located
at

    a_0=-b_0=±√(h̃_c/β̃)(2ϵ)^1/4.

The result of a numerical simulation of Eq. (<ref>) for
ϵ=0.005 is depicted in Fig. <ref>c,
showing how the displacement field inside the PTC begins to self organize
itself into the standing wave pattern predicted in Eqs. (<ref>)
and (<ref>). All simulations in this
paper were carried out using the Dedalus software package <cit.>.



 §.§ Emergent collective modes.


In the previous section, we demonstrated that the unstable, exponentially
growing modes of a nonlinear PTC herald the onset of a new phase.
After a short transitional period, the nonlinear PTC enters a steady
state with broken spatial translation and discrete time translation
symmetries. In this section, we focus on small fluctuation around
this new state. We demonstrate that the lattice structure of the electric
field supports soft, wavelike excitations, which propagate through
the lattice, similar to phonons in a crystal lattice. These excitations
are Goldstone-like modes that stem from the breaking of the continuous
spatial translation symmetry. Additionally, the steady state exhibits
a gapped mode, which is absent in the symmetric state.


These effects are most pronounced in the weakly damped regime,
where h≫ h_c, to which we will stick for the rest of this discussion.
We begin with the soft modes. Note that the phase of the spatial part
of the standing wave in Eq. (<ref>) is chosen spontaneously.
Replacing x→ x+ϕ/q^* in Eq. (<ref>),
we obtain another valid solution to Eq. (<ref>). Following
the usual logic of the Goldstone theorem <cit.>,
we expect that, since a uniform offset ϕ is inconsequential
for the system's dynamics, a long-wavelength spatial dependence of
ϕ will only have a small impact on the time evolution. It remains
to show that a long-wavelength perturbation of the form ϕ(x,t_0)=ϕ_Qe^iQx,
where Q≪ q^*, will propagate through the system as a wave with
frequency ω_G(Q), and that ω_G(Q)→0,
as Q→0, i.e. the dispersion ω_G(Q)
is soft. Indeed, we find that this is true (see Supplement): For |ϕ|≪1,
the propagation of the Goldstone-like mode is described by the equation

    1/c^2∂^2ϕ/∂ t^2=(1-1/3h)∂^2ϕ/∂ x^2.

The dispersion of the mode is ω_G(Q)=c_GQ,
where c_G=√(1-h/3)c. Fig. <ref>a
shows the numerical simulation of Eq. (<ref>), where a
boundary condition of the form ϕ(x=0,t)=-0.3sin(ω_G(Q)t)
was applied. The perturbation at the boundary propagates as a Goldstone-like
mode whose dynamics is described by Eq. (<ref>),
and forms a standing wave, such that the phase throughout the PTC
is given by ϕ(t,x)=[sin(ω_G(Q)t-Qx)+sin(ω_G(Q)t+Qx)]/2.
For the simulation, we chose Q=2π/Λ, where Λ=20λ
and λ=2π/q^*. We also demonstrate the propagation of
the Goldstone-like mode starting from an initial field D(x)=√(h̃/β̃)cos(q^*x-0.2sin(q^*x/5))
at t=0. The propagation of the distortion for long times is shown
in Fig. <ref>b.

In addition to the soft Goldstone mode, the steady-state of Eq. (<ref>)
supports spatially uniform oscillations of the amplitudes a and
b around the minima of Eq. (<ref>). This amplitude
oscillation can be viewed as a massive Higgs mode. To derive the oscillation
frequency, we expand the effective Hamiltonian (<ref>)
around these minima:

    H≈-h̃^2Ω/4β̃+h̃Ω(δ b^2+δ a^2).

Hamilton's equations read


    δḃ    =2h̃Ωδ a
    δȧ    =-2h̃Ωδ b.

Solving for a(t), we find δä+(2h̃Ω)^2δ a=0.
Thus, the PTC supports uniform amplitude oscillations of frequency
2h̃Ω. For finite γ, the above equations (<ref>)
will obtain a damping term and read δḃ=-γ b/2+2h̃Ωδ a,
δȧ=-γ a/2-2h̃Ωδ b. The solution
of these equations (including damping) is shown in Fig. <ref>b
together with the field amplitude D(t,2π/q^*). We
point out that if no dissipative term is included in Eq. (<ref>),
the PTC, starting with an initial field distribution, will not converge
to a fixed point. Instead, the amplitudes a, b will follow the
Hamiltonian dynamics of Eqs. (<ref>). 



 §.§ Symmetry breaking and pattern formation in two dimensions.



So far we have considered a one-dimensional PTC. A similar symmetry
breaking transition can be observed in higher dimensions. Here, the
spatial second derivative in Eq. (<ref>) is replaced by
a laplacian: ∂^2D/∂ x^2→∇^2D.
Focusing on the isotropic case, we omit vector signs. The analytical
calculation presented so far can be extended to higher dimensions.
There is, however, a crucial difference to consider. In one spatial
dimension, the only wavenumbers for which the resonance condition
(<ref>) is fulfilled are ± q^*. Plane
waves with these two wavenumbers arrange themselfs in a standing wave
pattern as shown in Fig <ref>c. In two dimensions,
all wavevectors lying on a circle of radius q^*=Ω/c around
the origin are resonant. From studies of Faraday waves, it is known
that parametrically excited waves often arrange themselves in stripes,
square patterns or hexagonal patterns. These arrangements correspond
to one, two or three critical wavevectors 𝐪^* that are
chosen spontaneously during the transition to the symmetry breaking
state <cit.>.

For the simple nonlinearity of Eq. (<ref>), we expect a
stripe-pattern <cit.> similar to the one-dimensional
scenario. On an infinite domain, the direction of 𝐪^*
will be chosen spontaneously. To show that the wavepattern for a two
dimensional PTC is indeed a stripe arrangement, we extend the ansatz
(<ref>) to N wavevectors 𝐪_i^*
with a relative angle θ_N, such that 𝐪_i·𝐪_i+1=𝐪_N·𝐪_1=q^*2cosθ_N
with θ_N=π/N. The driving term in Eq. (<ref>)
does not supply any momentum to the system. Therefore, each mode with
a wavevector 𝐪_i^* must be balanced by a mode with
-𝐪_i^*. Finally, we assume that all modes share the
same time dependent amplitude. Summing up, our ansatz is

    D=(a(t)cos(Ω t)+b(t)sin(Ω t))∑_i=1^Ncos(𝐪^*·𝐱).

When inserting this ansatz into Eq. (<ref>), the linear
terms can be treated similarly as in the 1D case. The nonlinear term,
however, requires additional care. As in the 1D case, we project the
nonlinear part onto the N resonant modes of the ansatz (<ref>)
and ignore fast oscillating, off-resonant terms. The equations of
motion for the amplitudes a and b are similar to Eqs. (<ref>)
where β has to be replaced by β̃_N=9β(2N-1)/32.
Consequently, we have to make the same replacement β̃→β̃_N
in the potential (<ref>). The new
potential has minima at a_N=-b_N=±√(h̃_c/β̃_N)(2ϵ)^1/4.
The depth of these minima depends on the number of modes N:

    V(b̅_N=a_N-b_N)∼1/β̃_N∼1/2N-1.

A similar conclusion can be reached in the weakly damped regime of
Eq. (<ref>). The depth of the minimia of the Hamiltonian
again behaves as ∼1/β̃_N. For more complex nonlinearlities
containing a richer structure of spatial derivatives, the dependence
of the nonlinear term on the angle θ can compensate for the
tendency of shallower minima for higher N. Typically θ=π/2,
or θ=π/3 is favored, leading to square, or hexagonal wave-lattices.
In general, one expects that the system converges to the deepest minimum
– in our case, to the N=1 minimum, leaving us with a single standing
wave. The orientation and phase are chosen spontaneously. 

Finally, we test our prediction of a stripe-pattern for the field
amplitude in the symmetry broken state by performing a simulation.
Starting from infinitesimal random noise on a square domain of size
10λ with periodic boundary conditions, we let the system
evolve in time and track the field amplitude. The results for the
near critical case with ϵ=0.3 are shown in Fig. <ref>.
Interestingly, before the PTC reaches the final stripe state, it goes
through a series of nearly-periodic patterns, which can be stable
for many (up to hundreds) oscillation periods.



§ DISCUSSION


Summing up, we investigated the long term behavior of a photonic time
crystal with a generic Kerr nonlinearity. We have shown that the initially
exponential growth of unstable gap modes in these systems heralds
the transition to a new symmetry breaking steady state. In this state,
the field configuration is a standing wave lattice, where the lattice
spacing is determined by the resonant wavenumber q^*, fulfilling
the resonance condition Ω=cq^*. The breaking of continuous
spatial translation symmetry is accompanied by the appearance of soft
Goldstone-like modes, which consist of propagating, wavelike lattice
distortions, similar to phonons in a crystal lattice. In addition,
gapped modes that correspond to amplitude oscillations of the electric
field emerge. The dispersions of these emergent modes depend on the
driving strength, which opens the possibility to use nonlinear PTCs
in the symmetry broken state as tunable metamaterials. We also investigated
the behavior of nonlinear PTCs in two spatial dimensions, and showed
that the electric fields self-organize to stripe patterns that break
the translational, and rotational symmetries of the system. This behavior
could help to create controllable, spatially structured high intensity
electric fields.

Furthermore, the predicted steady state exhibits subharmonic response,
and breaks the discrete time translation symmetry of the periodic
modulation. This is a common feature of many types of parametrically
driven nonlinear waves <cit.>,
as it has been pointed out in recent literature on time crystals <cit.>.
Thus the nonlinear PTC exhibits discrete time crystalline behavior
in the sense of a phase of matter <cit.>.

We thank S. Gom and M. Segev for useful discussions. E.K. acknowledges
financial support by the Helen Diller quantum center.









§ SUPPLEMENTARY MATERIAL
 












 §.§ Expansion close to the instability threshold. Continuity of the phase
transition.


Here, we explain the technical details involved in deriving the gradient
descent dynamics of Eq. (<ref>) in
the near-critical regime. First, we reformulate the amplitude equations
(<ref>) in terms of the small parameter ϵ=(h-h_c)/h_c=(h̃-h̃_c)/h̃_c:

    ȧ̅̇    =-Ωh̃_c(2+ϵ)a̅+β̃Ω/2b̅(a̅^2+b̅^2)
    ḃ̅̇    =Ωh̃_cϵb̅-β̃Ω/2a̅(a̅^2+b̅^2).

The small ϵ characterizes the growth of the unstable mode,
which defines a slow time scale. We introduce an expansion in terms
of ϵ, such that a̅=ϵ^1/4a̅_1+ϵ^3/4a̅_2...,
b̅=ϵ^1/4b̅_1+ϵ^3/4b̅_2+...
and ∂_t=∂_t_0+ϵ∂_t_1. At
orders ϵ^1/4 and ϵ^3/4, we obtain

    ∂a̅_1/∂ t_0    =-2Ωh̃_ca̅_1
    ∂a̅_2/∂ t_0    =-2Ωh̃_ca̅_2-β̃Ω/2b̅_1(a̅_1^2+b̅_1^2).

These equations describe a quick relaxation towards the fixed point

    a̅_1    =0
    a̅_2    =-β̃/4h̃_̃c̃b̅_1^3.

At order ϵ^5/4, we derive the equation

    ∂b̅_1/∂ t_1=h̃_cΩb̅_1+β̃Ω/2a̅_2b̅_1^2.

For a small ϵ, we can assume that the relaxation of a̅_2
towards the value given in Eq. (<ref>)
is infinitely fast, and write

    ∂b̅_1/∂ t_1=h̃_cΩb̅_1-β̃^2Ω/8h̃_̃c̃b̅_1^5.

Restoring the original variables t and b̅ we obtain

    ∂b̅/∂ t=ϵh̃_cΩb̅-β̃^2Ω/8h̃_̃c̃b̅^5.

This is the gradient descent dynamics described by Eq. (<ref>)
of the main text.



 §.§ Goldstone-like Modes


We add a time and space dependent phase to the steady state solution
found in the main text for low damping:


    D_ϕ(x,t)=±√(8h/9β)sin(Ω t)cos(q^*x+ϕ(x,t)).

Here, ϕ(x,t_0)≪1. We seek to derive an equation
governing the dynamics of the phase ϕ(x,t). To this
purpose, we write

    D_ϕ(x,t_0)     =A/2sin(Ω t-qx-ϕ(t,x))
         +A/2sin(Ω t+qx+ϕ(t,x)),

where A=√(8h/9β). In what follows, we will use the abbreviations
sin(Ω t-qx-ϕ(t,x))=sin(-)
and sin(Ω t+qx+ϕ(t,x))=sin(+)
and similar for cos-terms. Inserting Eq. (<ref>)
into Eq. (<ref>) and comparing the coefficients in front
of the cos(+) terms we find an equation for ϕ(t,x):

    1/c^2∂^2ϕ/∂ t^2=(1-1/3h)∂^2ϕ/∂ x^2.

The same equation follows from the cos(-) terms. The
sin(±) terms, on the other hand, give corrections
of order Q/q^* to the solution (<ref>).
These can be neglected for a smooth, small-Q perturbation of the
kind that we are considering. Compared to the velocity of electromagnetic
waves in the undriven medium, the velocity of the Goldstone mode c_G
is reduced by a factor that depends on the amplitude of the time varying
term in Eq. (<ref>):

    c_G=√(1-1/3h)c.





 §.§ Pattern formation in two dimensions




