Nonlinear two-level dynamics of quantum time crystals

A time crystal is a macroscopic quantum system in periodic motion in its ground state. In our experiments, two coupled time crystals consisting of spin-wave quasiparticles (magnons) form a macroscopic two-level system. The two levels evolve in time as determined intrinsically by a nonlinear feedback, allowing us to construct spontaneous two-level dynamics. In the course of a level crossing, magnons move from the ground level to the excited level driven by the Landau-Zener effect, combined with Rabi population oscillations. We demonstrate that magnon time crystals allow access to every aspect and detail of quantum-coherent interactions in a single run of the experiment. Our work opens an outlook for the detection of surface-bound Majorana fermions in the underlying superfluid system, and invites technological exploitation of coherent magnon phenomena – potentially even at room temperature.

P erpetual ground state motion in equilibrium defines a time crystal, but observing such motion is famously unfeasible 1 . Experimental time crystal realisations thus bend either the equilibrium [2][3][4] or the perpetuity 5-7 requirement, reaching stability only if isolated from the environment and the observer 1,[8][9][10] . Consequentially, coupling separate time crystals while retaining sufficient isolation is challenging, and time crystals have yet not been studied in a dynamic environment. We arrange spontaneous two-level dynamics of interacting time crystals, each consisting of 10 12 magnons, in the superfluid B phase of 3 He ( 3 He-B). In this system, the observable time crystal life time can be extended up to a thousand seconds 11 (10 9 periods of motion) in the absence of a driving force, while the underlying superfluid system provides intrinsic feedback for engineering coherent dynamics.
Magnons in 3 He-B arise as the quanta of transverse spin waves, associated with magnetisation that precesses about the external magnetic field H. At sufficient magnon density and low enough temperature, the precession synchronises spontaneously at uniform frequency ω and phase, forming a magnon Bose-Einstein condensate [12][13][14] . The spontaneous synchronisation can be demonstrated by pumping magnons to a higher energy level in the confining trap from which they spontaneously fall to the ground state 5,15 , or even by pumping incoherent magnons to the system using a noise drive 7,16 . This shows that the magnon state in the BEC is decoupled from the drive. The transverse spin precession of the magnon condensate, therefore, manifests the characteristic spontaneous periodic motion of a time crystal 5,6 .
The time crystal can be created using two different pumping techniques. Using a continuous drive yields a Floquet (discrete) time crystal. Here we use the pulsed technique where the drive is turned off before the time crystal evolution begins. This approach allows us to study uncontrived time crystal dynamics and interactions in the absence of external enforcement. The time crystal formation during the pumping pulse and its evolution thereafter is characterised by two timescales. The first timescale τ E~0 .1 s describes the time crystal thermalisation 17 , that is, how quickly the precession becomes coherent at the ground level in a trap, following the pumping of magnons. The second timescale τ N is the time crystal lifetime. In an isolated sample container τ N →∞ exponentially as temperature decreases. In practice there are also losses in the circuit that is coupled to the precessing spins for control and observation purposes. It is therefore necessary to allow for a finite τ N . The time crystal remains well defined as long as τ N ≫ τ E 5,6 .
In superfluid 3 He, Cooper pairs possess orbital momentum whose average distribution, parametrised by vector L, is axially symmetric in the sample container (Fig. 1). The time crystal is trapped in the middle of the superfluid sample by that distribution owing to spin-orbit interaction. We fine-tune the trap by adding a magnetic field profile as detailed in Methods [18][19][20] . We also place a free surface of the superfluid above the bulk trap centre 6 . The free surface distorts the distribution of L as shown in Fig. 2a, resulting in a second local minimum, located 3 mm above the bulk trap minimum. Magnons can be trapped and form time crystals in either of the traps or both of them simultaneously. In this Article we concentrate on the lowest energy level in each trap. We denote the time crystals "bulk" and "surface" corresponding to the physical location. The location of a time crystal is identified from experimental records by its response to changes in the magnetic field profile 6 .
Let us denote the bulk time crystal population (that is, the number of trapped magnons) N B and the surface population N S . The bulk and surface precession frequencies (ω B , ω S ) are determined by the profile of the confining trap, and the coupling Ω between the crystals by the overlap of their wave functions, as detailed in Methods. We will show that the dynamics of the coupled levels are described by the macroscopic two-level Hamiltonian where ℏ is the reduced Planck constant and t is time. Note that this two-level system is conveniently parametrised by a macroscopic Bloch sphere ( Fig. 2c) with the relative precession phase between the time crystals corresponding to the azimuthal angle and level populations corresponding to projection on the z axis, in a direct analogy with the Bloch sphere description of microscopic two-level systems such as qubits. In what follows we measure all frequencies in the frame rotating at the Larmor frequency ω 0 = |γH| taken at the centre of the bulk trap (γ ≈ −2⋅10 8 rad s −1 T −1 is the 3 He gyromagnetic ratio). The essential difference of Eq. (1) from a standard two-level Hamiltonian is the dependence ω B [N B (t)], which arises due to a nonlinear feedback provided by the spin-orbit trap: a large local magnon density expands the trap, changing the distribution of L, thus decreasing ω B . This mechanism is extensively studied in refs. 5,13,14,21 , and the outcome is schematically illustrated in Fig. 2b. Near the free surface the L distribution is fixed perpendicular to the surface. Therefore ω S is constant to a good approximation. As the time crystal populations decay, N B decreases and ω B increases. We can thus make the energy levels in the double trap cross by selecting suitable initial populations. A rigorous description of the macroscopic time crystal wave functions, the trapping potential, and the feedback mechanism are derived in Methods. We emphasise that all the relevant technical explanations can be found in Methods also where not explicitly referenced.
We note that another characteristic feature of time crystals, the lack of heating under continuous drive 22 , is also manifest in this system. Under continuous pumping, the number of magnons in Fig. 1 Schematic illustration of the experiment. The superfluid 3 He sample is contained in a quartz glass cylinder. The magnon time crystal (blue blob) is trapped in the middle of the container by the combined effect of a minimum in the static magnetic field, created using a pinch coil (green wire loop), and by the spatial distribution of the superfluid orbital momentum L (small green arrows). The coherent precession of magnetisation M (magenta cone) in the time crystal is observed using transverse pick-up coils. The static magnetic field H is oriented parallel to the axis of the cylinder. The ripple on the superfluid free surface is added for illustrational purposes. The two-level time crystal is schematically illustrated in Fig. 2. the time crystal is determined by the chemical potential that corresponds to the pumping frequency (see Methods). When this number would be exceeded, the time crystal spontaneously decouples from the drive 13 , thus preventing overheating. Similarily, during the slow population decrease after a pumping pulse, the chemical potential (precession frequency) is continuously adjusted to the changing magnon number 21 . Hence, the precession period and coherence become incommensurate with and thus independent of the drive pulse even if the drive was originally resonant.
In this Article, we study the dynamics of the two-level time crystal system, carrying out two experiments. In the first experiment, where the level crossing takes place at small N B , we showcase crossing dynamics that follow the textbook description: the system is initially in the ground state, but at an avoided level crossing both levels are populated owing to Landau-Zener population transfer. This "superposition" state continues to be modified by Rabi population oscillations after the crossing. The second experiment starts from a "superposition", the level crossing takes place at large N B where the feedback mechanism transforms into dynamically changing coupling Ω. Analysis of the second experiment shows that co-existing time crystals lay manybody interactions bare for the capable observer in a single run of the experiment. That is, in this case, the level crossing dynamics cannot be described analytically, but while coherent quantum phenomena are often hidden from direct inspection, time crystals have no such limitations.

Results
Basic two-level dynamics at small N B . The time crystal levels can be populated in desired proportion by a radio-frequency pulse via adjacent coils (Fig. 1). To highlight the two-level dynamics, we populate only the bulk time crystal in the beginning of the experiment shown in Fig. 3a. After the pulse, the coherent precession of magnetisation induces an oscillating signal in the coils, which allows inferring the precession frequency, and the signal amplitude yields the magnon number. These quantities are extracted from the experiment in Fig. 4. The pumping is followed by exponential decay of N B ðtÞ ¼ N B ðt ¼ 0Þ expðÀt=τ B Þ with time constant τ B , controlled by temperature as detailed in Methods.
In Fig. 3, the ground level is initially located in the bulk trap, N B decays at a rate determined by τ B , and ω B increases slowly as the trap recovers a narrower shape. Meanwhile, ω S remains constant. Hence, ω B eventually crosses ω S before levelling out. In a coupled two-level system, a level crossing has specific consequences: The observed frequencies are the (dressed) eigenfrequencies of the Hamiltonian which deviate from the undressed frequencies ω B , ω S in the Rabi regime Ω > |ω B − ω S |. Due to this hybridisation, the observed levels avoid crossing each other, and the global ground level smoothly switches from bulk to a b c Fig. 2 Time crystal two-level system. a The distribution of L (green arrows) confines magnons in two local minima, hosting two adjacent time crystals: one in the bulk of the superfluid (blue blob) and the other one touching the free surface (red blob). In each time crystal, magnetisation is precessing coherently, which couples to measurement circuitry as shown in Fig. 1. b Magnons in the bulk modify the confining trap created by the L distribution. When the bulk population is large (cyan blob), the textural trap is widened (red arrows), which modifies also surface time crystal's wave function (magenta blob). This increases the coupling between the states. Changes in the trap and the wave functions have been exaggerated for illustrational purposes. c The state of the two-level system (red arrow) can be illustrated using a Bloch sphere where radial distance corresponds to magnon number N B + N S , the relative phase between the time crystals' precession corresponds to the azimuthal angle ϕ, and the polar angle θ describes the relative weights of the two-level basis states in the "superposition" (see Methods). Frequency is plotted in the rotating frame Δω = ω − ω 0 . The excitation pulse at t = 0 is framed out for clarity.
Initially ω B < ω S , but as the population in the bulk trap decays, at t = 3.3 s, the global ground state moves to the surface in an avoided crossing. The excited state, now located in the bulk, is simultaneously populated. Rabi (Josephson) population oscillations are seen as a side band. Coupling extracted from the side band extrapolates to Ω/(2π) = 1.7 ± 0.4 Hz at the crossing, in good agreement with the fitted simulation value Ω/(2π) ≈ 1.4 Hz. b The numerical simulation recreates the population transfer and the side band, confirming that the population transfer is due to Landau-Zener transition (analysis in Fig. 4). In the absence of measurement noise also the side band of the surface time trace is weakly visible. c Subtracting the simulation from the experiment shows that the point-wise residuals remain smaller than 5%. The relative difference is normalised by the total signal at the crossing. In this measurement temperature was 180 μK and ω 0 /(2π) = 833 kHz. surface (from ω B to ω S ) as seen in Fig. 3a. Population transfer between the levels is also observed, as both of the levels are populated after the avoided crossing.
Traversing the avoided crossing adiabatically would allow the entire magnon population to follow the global ground state, but here some magnons move to the state with higher eigenenergy (precession frequency). This process is generally known as Landau-Zener-Stueckelberg-Majorana tunnelling, or Landau-Zener tunnelling. The population transferred depends on the rate of level crossing d|ω S − ω B |/dt at ω B = ω S . For usual coupled non-linear oscillators with damping, the crossing rate is determined directly by the damping. Here the decay rate τ B = 3.5 s corresponds to d|ω B [N B (t)]/(2π)|/dt2 4 Hz s −1 . Using this and the coupling directly extracted from the experiment as explained below yields the predicted Landau-Zener population transfer fraction 0.9% to the excited state, which is two orders of magnitude smaller than that observed in the experiment. In this case, only one of the two levels would be visible in Fig. 3a after the crossing. In similar experimental runs with smaller d|ω S − ω B |/dt we have observed population transfer up to 20 orders of magnitude larger than the corresponding Landau-Zener prediction.
We can analyse this striking mismatch by simulating the time evolution of the two-level Hamiltonian numerically. We feed the experimentally-determined bulk and surface decay rates, τ B and τ S , the corresponding initial populations, and the measured ω B [N B ] dependence to a numerical simulation of the two-level Hamiltonian (see Methods). The coupling constant Ω is used as a fitting parameter yielding Ω/(2π) = 1.4 Hz.
The outcome of the simulation, plotted in the same way as the experimental signal, is shown in Fig. 3b. We can directly compare it with the experiment by subtracting the simulated time-dependent Fourier spectrum from the experimental one (Fig. 3c). The simulation underestimates the change of the surface time crystal frequency near the crossing as shown in Fig. 4a,b, which causes the largest deviation between the Fourier spectra. Otherwise the typical deviation between the two signals, as normalised by the total signal at the crossing, is less than 5%. In particular, the simulation replicates the magnitude of the population transfer, that is, 60% of magnons move to the excited state. The simulated population dynamics are compared directly with the measured signal in Fig. 4c. We emphasise that repeated runs of the simulation with perturbed input parameters reveal that this qualitative level of population transfer is insensitive to the precise value of any of the input parameters.
To explain this observation, we extract the undressed frequencies from the simulation in the region near the avoided crossing (Fig. 4b). As a result of the feedback ω B (N B ), the bulk frequency is changing both owing to the slow decay of N B and because of the population transfer from N B to N S by Rabi oscillations. Thus, their combined effect increases the crossing rate d|ω S − ω B |/dt. The magnitude of the Landau-Zener population transfer is determined within~50 ms of the crossing 23 , and in this window ω B (t) can be linearised. Inserting the accelerated crossing rate to the Landau-Zener formula yields the expected population transfer of 61%, in good agreement with the simulated population transfer, 60% (Fig. 4d). That is, the population transfer follows the Landau-Zener description with the crossing rate taken at the instant of level crossing. Note that the Landau-Zener description is thus valid even if the crossing rate is regulated by intrinsic feedback. We conclude that the observed population transfer strongly supports the two-level interpretation of the time-crystal dynamics. Far from the avoided crossing, the two-level interaction is characterised by AC Josephson population oscillations between the levels 6 . Owing to the feedback in the bulk trap, the oscillations result in a side band that follows the bulk trace. The frequency of the population oscillations is set by the difference of the time crystal precession frequencies, equal to the difference of their chemical potentials 6 . Thus, the side band is separated from the bulk trace by |ω B −ω S |, as seen in Fig. 3a (see derivation in Methods). The amplitude of the population oscillations is determined by the coupling Ω, and the relative side band amplitude by the slope of ω B (N B ) (formulas are given in Methods). Thus, we can extract the coupling directly from the experimental data, yielding Ω/ (2π) = 1.7 ± 0.4 Hz in the level crossing region, in good agreement with the simulation fitted value Ω/(2π) = 1.4 Hz.
Put together, the Josephon population oscillations, the Landau-Zener population transfer, and the agreement on the two-level coupling independently extracted from the different aspects of the population dynamics confirm that the two-time crystals form a macroscopic two-level system.
Dynamical coupling regime at large N B . Magnon time crystal dynamics, enhanced by the nonlinear feedback, can be analysed also directly without resorting to a numerical simulation of the system. This is advantageous as it will allow untangling interactions involving multiple time crystals that go beyond the twolevel description. As a simple demonstration of this capability, we introduce a level crossing in a region where N B is an order of magnitude larger than above. The resulting trap modification affects not only ω B (N B ) but also the constriction between the time crystals, as sketched in Fig. 2b. This causes the coupling Ω to change dynamically in the course of the crossing. Both levels are populated in the beginning of the experiment (Fig. 5a) to allow following their dynamics directly.
In this experiment the coupling is changing, but qualitatively the dynamics follow a similar pattern as above: The ground state moves from bulk to surface when the undressed frequencies cross at t ≈ 3.8 s (Fig. 5b). The moment of the crossing is identified by a sharp increase in the ground state trace relaxation rate from τ À1 B ¼ 0:06 s À1 to τ À1 S ¼ 0:53 s À1 (Fig. 5c). The increase is attributed to increased dissipation in the surface trap due to surface-mediated emission of other spin wave modes 24 and potentially surface-bound Majorana states 25,26 , but a detailed study is left for the future. We note that the states can be identified also by adjusting the magnetic field profile and rerunning the experiment; the relaxation rate is a convenient shortcut for distinguishing the two levels.
Josephson population oscillations between the two time crystals are seen as the side band that follows the bulk time crystal trace. As explained above, the side band is separated from the bulk trace by |ω B − ω S |. This separation is characteristic of the Josephson effect, and it changes in time because ω B changes, as shown in Fig. 5d. The surface time crystal is not followed by a similar side band, because the surface trap is rigid and hence population oscillations result in no side bands (see Methods). A second bulk trace side band should be located symmetrically on the other side of the bulk trace, but it exactly coincides with the surface trace and is therefore not resolvable.
The side band amplitude allows us to extract the coupling between the time crystals (Fig. 5e). The extracted coupling is the largest in the beginning of the experiment and decreases when N B decreases. That is, the constriction between the time crystals is affected by the bulk trap modification (Fig. 2b), which makes the coupling larger when N B is large. This is qualitatively in line with the trap modification mechanism discussed in refs. 13,14,21 .
Near the avoided crossing interference effects prohibit direct access to the population oscillation in the experiment. However, in a two-level system the coupling can also be extracted from the minimum frequency separation of the dressed frequencies of the two levels at the avoided crossing, which is equal to 2Ω. This is done by interpolation in Fig. 5a. The result is shown by the horizontal line in Fig. 5e, in good agreement with the dependence extracted from the side bands. Note that near the avoided crossing the separation of the dressed (observed) frequencies equal to the difference between the Josephson frequency and the Rabi frequency. That is, Rabi population oscillations smoothly replace the Josephson oscillations, increasing the population oscillation frequency as compared with the Josephson frequency, which goes to zero.
It is worth noting that the relaxation rate of the bulk time crystal depends on whether it is the global ground state or the global excited state. In Fig. 5c the bulk relaxation time is τ À1 B ¼ 0:06 s À1 until the level crossing (ground state), and it increases to τ À1 B ¼ 0:2 s À1 after the level crossing (excited state). The same observation, correspondingly, applies to the surface time crystal relaxation rates. We emphasise that such change is never observed in the absence of the level crossing, for example, if the bulk time crystal is the ground state throughout the experiment. That is, the excited state seems to slowly leak magnons to the ground state. This observation hints that there is an additional incoherent channel that allows magnons to move from the excited state to the ground state, showing independently that the two time crystals interact and that the level crossing has physical consequences that penetrate the dynamics in the two-level system.
The above analysis confirms that the two-level description is valid and robust against dynamic variation its parameters, and that direct experimental observations provide continuous access to all relevant aspects of the interaction.

Discussion
To summarise, we have shown that the dynamics and interactions of the two adjacent magnon time crystals are quantitatively described by a two-level Hamiltonian. The levels are modified by a nonlinear feedback, arising owing to spin-orbit interaction in the underlying superfluid system. This allows engineering intrinsic time crystal dynamics in the absence of continuous external drive. We show that when the two-level eigenfrequencies approach one another, the coupling between the levels results in an avoided crossing with ensuing Landau-Zener population transfer from the global ground state to the excited state. Rabi population oscillations, combined with the feedback mechanism, increase the population transfer by orders of magnitude. This is quantified by comparing numerically simulated population dynamics with the experiments. We also show that all relevant observables and parameters including the eigenfrequencies and the coupling between the time crystals can be simultaneously extracted from the experiment. We emphasise that each measurement sequence shown in this Article corresponds to a single run of the experiment, but the phenomena are well reproducible.
We have shown that the spin-orbit interaction can be harnessed to create a nonlinear feedback for magnons in a coherent time crystal system. Nonlinear feedback is needed for spin-based versions of quintessential quantum devices such as the SQUID. It remains an interesting task to explore the time-crystal two-level system further by demonstrating parametric pumping of magnons and logic gate operations between the two levels. For example, parametric pumping can be arranged by modulating the magnetic part of the trap at frequency Ω. Additionally, any number of co-existing time crystals can be accommodated in a magnetic landscape to increase the number of degrees of freedom, and the flexible trap can be turned off by adjusting the external magnetic field. These are important capabilities for realising magnon-based devices [27][28][29][30][31][32] . To access phenomena such as quantum entanglement, few-magnon operations can be implemented using nano-fluidic confinement and ultra-sensitive NMR techniques 33,34 . We emphasise that similar physical phenomena including quasiparticle Bose-Einstein condensation and the emergence of time crystals can be accessed in certain solid-state room-temperature systems, for example based on magnons in YIG films [35][36][37][38][39][40][41][42] . This opens the outlook of quasiparticle-based coherent on-chip applications in ambient conditions, including coherent quantum information processing 27,28,[30][31][32]39 .
The three-dimensional topological superfluid is wrapped by a two-dimensional system of surface-bound quasiparticles, among them Majorana fermions [43][44][45][46][47][48][49] . At the free surface there are no impurities (unlike at sample container walls), and surface-bound Majorana fermions are expected to manifest themselves as detectable zero-temperature magnetic dissipation 25,26 . Majorana fermions have remained elusive despite a decade of searching in different condensed matter systems 50 . The hybridised two-level state is in direct contact with the free surface of the superfluid, making and extremely sensitive probe for the bound Majorana states. We have provided preliminary evidence that the surface induces magnetic dissipation and that details of this signature can be explored using the time-crystal two-level system.

Methods
Experiment. The superfluid 3 He sample is placed in a cylindrical quartz-glass container (15 cm long, 6 mm diameter) in a nuclear demagnetisation refrigerator (Fig. 1). The lower end of the sample container connects to a volume of sintered silver powder surfaces, thermally linked to the nuclear refrigerant. This allows cooling the 3 He down to 130 μK. Temperature of the superfluid is measured using a quartz tuning fork 51,52 , and pressure is equal to saturated vapour pressure, which is vanishingly small at these low temperatures. The superfluid transition temperature at saturated vapour pressure is T c ≈ 0.9 mK. The sample container is surrounded by two transverse NMR coils, which are part of a tank circuit resonator with Q ≈ 150, and a pinch coil used to create an axial minimum of the magnetic field. The resonance frequency of the tank circuit can be tuned in eight equidistant steps between 550 kHz and 833 kHz, corresponding to external magnetic fields between 16.5 mT and 25 mT. The signal is amplified by a cold preamplifier 53 and room temperature amplifiers.
The free surface is located 3 mm above the centre of the magnetic field minimum. The location of the free surface is adjusted by removing 3 He slowly until the desired location is achieved, measuring the pressure of 3 He gas in a calibrated volume that results from the removal of liquid from the originally fully filled sample container. The outcome is favourably compared with the observed magnon spectrum and a numerical model of the trap. The resulting two traps for magnons are detailed in the next Section.
The time crystal wave function can be written as Ψ = ae −iωt , where t is time, ω is the precession frequency related to the chemical potential μ = ℏω, the phase term e iφ is contained in a, and the number of magnons N = |a| 2 . The tipping angle of the precessing magnetisation β M , measured from the magnetic field H, parametrises the spatial profile of the wave function, N ¼ jaj 2 / R sin 2 β M 2 dV. The signal induced in the pick-up coils (Fig. 1) is sinusoidal, corresponding to the magnetisation along the axis of the NMR coil, or in other words, real part of the rotating complex wave function, e −iωt . The measured signal amplitude is proportional to the amplitude of the time crystal wave function, where c contains the so-called filling factor of the state within the NMR coils, the amplification provided by the tank circuit resonator and other amplifiers in the measurement circuit 53 , and physical constants 20,21 . A desired level in the trap can be populated by a radio-frequency pulse via the pickup coils, followed by slow population decay owing to two mechanisms: The fermionic thermal excitations of the superfluid cause non-hydrodynamic spin diffusion 18 . This contribution can be made exponentially small in the zero-temperature limit (1000 s life time has been achieved 11 ) or dominant at higher temperatures. Observing and controlling the quasi-perpetual time crystal motion inevitably causes also external dissipation, 54 in our case radiation losses in the measurement circuitry 18 . Both of these dissipation mechanisms cause exponential population decay in time, in combination described by time constant τ N . The time crystals are well defined provided the life time, which here is τ N~1 0 s, is much longer than the time it takes for the time crystal to form after the pulse (here τ E~0 .1 s) 5,6 . The two level system, in the absence of coupling between the states, is described by the "superposition" wave function Only the relative phase enters the dynamics of the system. Hence, b can be chosen to be real, and the combination b, s conveniently illustrated by a macroscopic Bloch sphere (Fig. 2c): The surface corresponds to states with total magnon number N 0 = |b| 2 + |s| 2 = N B + N S , and the interior to smaller magnon numbers reached during the population decay. The weights of the basis states in the superposition, that is, the fraction of the total population in the bulk (surface) state is given by the polar angle θ with N B ¼ N 0 cosðθ=2Þ (N S ¼ N 0 sinðθ=2Þ). The relative phase ϕ corresponds to the azimuthal angle in the x-y plane of the sphere. It evolves in time according to We note that controlling the relative phase is beyond the scope of the present work and requires adjusting the NMR coil geometry.
Level dynamics in a flexible trap. 3 He-B is a p-wave superfluid, hence, the orbital momentum of the Cooper pairs is equal to one. In the sample container cylinder, the average orbital momentum L is distributed symmetrically ("texture", Fig. 1) owing to the orienting effects of the magnetic field and the container walls. In addition, we create an axial minimum of H using a pinch coil, which confines the magnons due to the Zeeman energy. The bulk trapping potential U(r) = U H + U L therefore has a magnetic part, and a component created by the L distribution owing to the spin-orbit interaction Here ω 0 (r) = |γH(r)| is the local Larmor frequency which depends on position r, Ω B is the B-phase Leggett frequency, γ ≈ −2⋅10 8 rad s −1 T −1 is the gyromagnetic ratio of 3 He, and the order parameter distribution is parametrised by the tipping angle of the orbital anisotropy axis, β L (r), measured from the direction of the magnetic field H, oriented along the cylinder axis.
Bringing the free surface above the trap centre distorts the order-parameter trap as β L = 0 at the free surface, creating a local minimum at the surface. Note that we study the time crystals in a frame rotating at the Larmor frequency ω 0 where the uniform magnetic field is absent. Where the notation ω 0 is used without an explicit reference to position, this means Larmor frequency in the middle of the bulk trap, corresponding to the minimum of the harmonic trapping potential. The time crystals located in the two traps can be identified and their frequencies adjusted by changing the profile of the field minimum, aided by the different relaxation rates. Below we concentrate on studying the feedback created by the flexible bulk trap.
The harmonic bulk trap has a radial trapping frequency ω r /(2π)~200 Hz corresponding to U L and an axial trapping frequency ω z /(2π)~20 Hz corresponding to U H . The resulting precession frequency is ω 0 + ω r + ω z /2. Therefore the axial trap can be neglected in the below analysis. It is thus convenient to measure all frequencies in the frame rotating at ω 0 . A more detailed analysis of the bulk trap can be found in refs. [18][19][20] .
The textural part of the trapping potential feels local magnon density due to spin-orbit interaction: The equilibrium texture minimises a range of free-energy contributions, including the orienting effects of the magnetic field and the sample container walls 55 . An important additional contribution is the spin-orbit interaction energy where ΨðrÞ / sin 2 β M =2 contains the spatial variation of magnon density which gives rise to the feedback effect. That is, the bulk trap profile and the shape of the time crystal wave function depend on N B so that dω B /dN B < 0. In the limit of large magnon number the bulk trapping frequency follows 13 Here k > 0 depends on the rigidity of the textural trap and the profile of the magnetic field minimum, p ≈ 5/7, 13 and ω B stands for the time crystal trapping frequency in the limit of zero magnons. We emphasise that although ω B changes during the decay of the magnon time crystal, the change is very slow as compared with ω 0 /(2π)~1MHz, and we can thus assume that the wave functions always correspond to the instantaneous trap shape 21 . Note that the surface trap is rigidified by the adjacency of the free surface, and ω S is therefore independent of N S to a good approximation.
It is possible to describe the self-trapping effect numerically in a self-consistent calculation of the order parameter texture 55,56 , the resulting trap 20 , the time crystal wave function 13,14,21,57 , and population decay 18,19,24,58 . That is however not necessary for understanding the experiments presented in this Article, because finding a general form of Eq. (7) can be circumvented by fitting and numerical differentiation of the experimental data where necessary, and all other effects can be measured independently. For simplicity, we refer to Eq. (7) in the below discussion, but the reader should bear in mind that the general form of the nonlinearity is more complicated.
Josephson coupling. Let us study the observable consequences of the population oscillation. We use the language of the Josephson effect, analogous to the AC Josephson effect 6 , as the oscillation amplitude can only be reliably extracted from the experiment far from the avoided crossing. Near the avoided crossing one should use the more general Rabi oscillation picture.
The amplitude of the AC Josephson population oscillation is Here Ω is the coupling, and ΔN B ≡ −ΔN S . The Josephson frequency is ω J = |ω B − ω S |. This oscillation modulates the bulk condensate frequency ω B as follows from the self-trapping Eq. (7). The frequency modulation (FM) is sinusoidal to a good approximation. This is because the amplitude of the population oscillation is small as compared with the total population, and Eq. (7) can be linearised.
The resulting instantaneous bulk time crystal frequencyω B can be written as Here Δω B is the FM amplitude. It is connected to the population oscillation amplitude ΔN B by Fourier decomposition of the resulting frequency-modulated signal yields Here J n is the Bessel function of the first kind of order n. The bulk main trace corresponds to n = 0. Combining the above expressions, and denoting the first side band (|n| = 1) amplitude as A SB , the coupling term can be linearised and expressed in quantities that can be directly measured: Here we assumed that the filling factors of the bulk and surface states in Eq. (2) are equal and constant. Where the time crystal shapes are changing due to changes in the trap profile, the coupling extracted using the above expression is therefore only approximate.
The side band of the bulk time crystal is seen in Fourier analysis of the experimental signal (Fig. 5d). The coupling extracted from this record using Eq. (12) extrapolates to Ω/(2π) ≈ 1.7 Hz at the crossing, in good agreement with the fitted simulation value Ω/(2π) ≈ 1.4 Hz. Note that there should be another side band symmetrically at lower frequency than the bulk trace, but it is covered by the surface trace at exactly the same frequency.
The surface trap is only weakly modified in similar fashion, yielding no visible side bands in the experiment. That is, the AC Josephson effect in a fully rigid trap results in no side bands owing to complex interference of the two wave functions. This can be confirmed by solving the dynamics of the rigid non-decaying coupled system analytically. We used this result to test the validity of the numerical NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-30783-w ARTICLE NATURE COMMUNICATIONS | (2022) 13:3090 | https://doi.org/10.1038/s41467-022-30783-w | www.nature.com/naturecommunications simulation discussed below. Note that Ref. 6 misleadingly implies that population oscillations directly cause the side bands even in the absence of nonlinear feedback.
Near the avoided crossing one should use the more general Rabi oscillation picture. Solving for the eigenfrequencies of the Hamiltonian in the Rabi regime yields the Rabi frequency ω R ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi ðω B À ω S Þ 2 þ ð2ΩÞ 2 q . In the limit Ω ≪ |ω B − ω S | this is reduced to ω R = ω J . The region where ω R ≠ ω J is not directly visible in the experiments due to interference effects.
Landau-Zener tunnelling. In the presence of exponential population dissipation, the time crystal population follows where 1/τ α is the relaxation rate of the measured signal (2), and α is either B for the bulk or S for surface. For the surface time crystal this makes little difference other than that the population decays slowly. The bulk time crystal frequency ω B depends on N B according to Eq. (7), and the frequency therefore increases during the decay. Hence, we have obtained the flexible two-level system described by the Hamiltonian (1). Let us choose ω B (N B = 0) > ω S and N B (t = 0) such that ω B (N B (t = 0)) < ω S . Now the frequencies of the surface and bulk time crystals will cross in the eigenbasis where Ω = 0. If Ω > 0 and N B decreases adiabatically, magnons in the bulk trap will smoothly move to the surface trap, remaining in the global ground state in an avoided crossing. The minimum frequency separation of the global ground state and the excited state at the avoided crossing is 2 Ω, as can be solved from the Hamiltonian.
If the avoided crossing is passed non-adiabatically, a part of the ground state magnon population moves to the excited state. This phenomenon is known as the Landau-Zener-Stueckelberg-Majorana effect. In our case this means that after the avoided crossing some population remains in the bulk trap, which corresponds to the new excited state in the system. The fraction of population promoted to the excited state is 59 where ∂ t stands for time derivative. Note that while in the canonical Landau-Zener problem the time derivative is constant, in our case it keeps changing. However, the magnitude of the Landau-Zener population transfer is determined within a time window $ 1= ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi j∂ t ðω B À ω S Þj p of the level crossing 23 (≲100 ms wide in our experiment). Therefore, ω B (t) can be linearised, or in other words, the time derivative taken at the avoided crossing gives the correct Landau-Zener population transfer.
Numerical simulation. The time crystal two-level Hamiltonian can be combined with the slow decay into a pair of equations: Here the right hand side corresponds to the Hamiltonian (1), and i is the complex unit. This pair of equations can be solved numerically. Our main motivation for the numerical simulation is to show that the simple two-level Hamiltonian describes the dynamics of the system exhaustively, that is, that the large population transfer is explained by the intrinsic dynamics of the two-level Hamiltonian. The most important test for this picture is the avoided crossing and the related population transfer. Reproducing the Rabi oscillations that result in a frequency-modulation of the time crystal frequencies is a secondary test.
The initial time crystal wave functions, the time crystal decay rates, and the bulk trap self-trapping power law, Eq. (7), can be extracted from experimental data independently and used as the parameters of the numerical simulation. The coupling Ω at the avoided crossing cannot be directly extracted from the experiment, and is used as a fitting parameter. To compare with the measured signal we also need the filling factors c α . They are used as fitting parameters as well.
We note that to reproduce the experimental signals in general, three additional effects need to be included in the simulation: (i) The surface time crystal frequency depends on the population in the bulk trap and (ii) also on the population in the surface trap; ω S = ω S (N B , N S ). (iii) Both τ B and τ S change at the avoided crossing, and to allow for this transition to take place smoothly, we use a smooth interpolation function between the asymptotic relaxation rates in the Rabi regime with the width of the crossing region being another fitting parameter with valuẽ Ω −1 . Effect (i) is due to the widening of the constriction that separates the time crystals (Fig. 2). This connection is included in the simulation, and it is seen as a decrease of ω S in Fig. 4a at t < 1 s, where bulk population is large and surface population negligible. However, this effect can be safely neglected in the analysis of the Landau-Zener effect in Fig. 3, because the avoided crossing takes place at small N B . The second dependence (ii) is what produces the frequency-oscillations of the magenta line in Fig. 4c. Both effects can be extracted from experimental data independently.

Data availability
The data used in this study are available in the Zenodo database under accession code https://doi.org/10.5281/zenodo.6510863.

Code availability
The simulation codes and guidance in their use are available from the corresponding author upon reasonable request.