Breathing mode dynamics of coupled three-dimensional chiral bobbers

Magnetic skyrmions are topological, nanoscale spin structures, which are embedded in a ferromagnetic background, with their cores being polarized oppositely to the background. In general, they can be categorized by the symmetry of their domain walls into Néel- or Bloch-type skyrmions. Their chirality is determined by the Dzyaloshinskii–Moriya interaction (DMI), which can be present in non-centrosymmetric systems and which will result in left- or right-handed skyrmions. The specific symmetry of the DMI is related to its origin via the type of symmetry breaking within the system. Owing to their topological nature, skyrmions are particularly stable. Ever since skyrmions have been brought into focus as spintronic information carriers, a plethora of scientific activities have shed light on their static and dynamic behaviors.^1–4 Experimentally, Bloch skyrmions have been observed in a wide range of B20 materials,^5–7 while Néel skyrmions were predominantly observed and studied in ultrathin multilayers with broken inversion symmetry at the interfaces.^8–13 Note that Néel skyrmions were also reported in bulk systems with specific crystal symmetries^14,15 and as surface states of B20 crystals.¹⁶ 

Parallel to experimental research, a zoo of topological objects have been theoretically predicted, some of which emerge when extending the planar description of a skyrmionic excitation to the third dimension.¹⁷ The simple stacking of skyrmions of the same size and chirality will, e.g., generate a skyrmion tube (also called a string).^18–21 More complex structures can be produced by twisting the string-like skyrmion or by connecting the ends to form a closed loop, which are often referred to as hopfions.^22–26 Such knotted magnetic structures can also be formed in systems with higher order magnetic interactions,²⁷ as well as in models of frustrated magnets, where they have been first predicted.²⁸ 

Another three-dimensional topological soliton is the chiral bobber.^29–32 It can be imagined as a stack of skyrmions that continuously decrease in size, eventually terminating in a Bloch Point (BP)—a magnetic point singularity that can also, intriguingly, be described as a point-like quasi-monopole. In theory, bobbers can be readily generated by interfacing skyrmion-hosting materials with a thin magnetic layer with high perpendicular magnetic anisotropy and high exchange stiffness. The interaction of this layer with the skyrmion hosting material causes a pinching off of the skyrmion tube. Most commonly, the additional layer is simply treated as a pinned layer with fixed uniform magnetization.^22,24 When a homogeneous skyrmion string is pinched off at both ends, two terminating BPs can be generated. These structures are referred to as torons,^24,33 monopole strings,²² or dipole strings.^34,35 Note that the term dipole string is only used when referring to an object with two BPs of opposite topological charge. Another way to interpret this structure is by imagining a stack of two bobbers. In our work, this picture will be the most appropriate due to its capacity to describe the ability of two bobbers to interact. We will thus refer to these structures as bi-bobbers.

Experimental evidence by real-space imaging of bobbers was reported by the means of off-axis electron holography by comparison of skyrmion and bobber contrast.³⁰ Reciprocal-space imaging to distinguish tube lattices penetrating through the entire specimen and bobber lattices with a termination located well inside the material is reported by quasi-tomographic, resonant, elastic x-ray scattering.³² Furthermore, a distinction by topological Hall effect measurements is also suggested.³¹ Nevertheless, the unambiguous experimental distinction between skyrmion tubes and bobbers or bi-bobbers in systems that cannot be substantially modified to be suitable for resolving the nanoscale magnetization in advanced imaging techniques (e.g., in TEM) is rather challenging. Additionally, in case of isolated objects or missing references for the contrast of bobber and tube, the identification may also be difficult. Most commonly available techniques utilizing static imaging still cannot unambiguously distinguish between the distinct topological fingerprints of more complex objects, as they mostly differ in their projected 2D images, i.e., without the information about the third dimension, so that many topological objects can lead to the same 2D contrast. Especially, spin textures with differing in-plane components will show the same result with 2D projected imaging techniques, but might have different internal modes, as will be useful in another method to distinguish spin textures explained in the next paragraph. It should be noted that with the prospect of new 3D imaging techniques, such as laminography,^36,37 these issues may be overcome in the near future, but this would require large-scale facilities.

On the other hand, a novel and promising all-electrical dynamic method allows for the detection of different magnetic structures based on their unique internal resonance modes.³⁸ These modes can be excited by periodic external stimuli in the microwave range via a magnetic field,^22,24,39,40 spin torques,^40,41 light,⁴² or even by moving skyrmions in nanotracks with periodic constrictions.⁴³ Internal modes of skyrmionic excitations can be classified as breathing or gyration modes. The former, which will be relevant for this work, is a periodic change of the structure’s size resembling the motion of breathing lungs. Depending on the specific magnetic structure, the internal modes can differ in the number of existing resonances, their resonant frequencies, and, most importantly, the dispersion relations as a function of material parameters or, more pertinently, as a function of external fields. Experimental studies of internal modes of skyrmion strings have been successfully conducted using coplanar waveguides and network analyzer ferromagnetic resonance (FMR) setups.⁴⁴ Micromagnetic simulation studies on the breathing of skyrmions,⁴⁵ antiskyrmions,⁴⁶ and synthetic antiferromagnet skyrmions⁴⁰ have been reported for frequencies up to 100 GHz. Recently, toron and hopfion breathing modes have been investigated up to 15 GHz.²⁴ 

In this work, a new type of bi-bobber living in a synthetic antiferromagnet (SAF) is introduced. These new bobber systems thus yield additional degrees of freedom in their stacking, which we study with focus on the effects of different chiralities, finding a dependence of the skyrmion stray field on its winding direction. Such a dependence of the strength of a magnetic stray field on the orientation of the magnetic elements producing it is known from certain arrangements of permanent magnets, usually referred to as the so-called Halbach effect.^47,48 It has profound effects on the properties of bi-bobbers. Furthermore, we will compare breathing modes of bi-bobbers and skyrmions in SAFs, excited by periodic out-of-plane magnetic fields with frequencies up to 60 GHz.

Isolated skyrmionic excitations can be readily studied in disks due to the inherent symmetry of edge states and stray field effects. Thus, we have chosen disk geometry with a diameter of 100 nm and height of 20 nm, which can comfortably accommodate our magnetic textures. Using a cell-size of 1.5625 × 1.5625 × 1 nm³, a total of 64 × 64 × 20 cells were simulated. For the chiral ferromagnet, an exchange stiffness of A_ex = 7 pJ m⁻¹, a saturation magnetization of M_s = 450 kA m⁻¹, and an interfacial DMI of D = 2.65 mJ m⁻² were chosen, lying between values of other reports for the purpose of comparability.^24,40,45 A perpendicular magnetization anisotropy (PMA) was achieved by setting the out-of-plane uniaxial anisotropy constant to K_eff = 50 kJ m⁻³. To stabilize small, geometrically constrained bi-bobbers, we added uniformly magnetized layers at the top and bottom of the simulated stack. For this work, we decided not to fix these layers artificially but, instead, added five layers to the top and bottom with stronger PMA (K_eff = 1 MJ m⁻³) and a slightly higher exchange stiffness (A_ex = 20 pJ m⁻¹). For simulations of the skyrmion tubes, the uniform top and bottom layers were not included to encourage homogeneous tube formation. Also, for better comparability, exchange stiffness and PMA were slightly increased when simulating skyrmion tubes to achieve similar sized tube diameters compared to bobber diameters at the interface.⁴⁹ The Gilbert damping coefficient was set as α = 0.01 to enable the response in the microwave frequency range.

With this configuration, especially by utilizing an interfacial DMI (iDMI), the simulated model can be thought of as thin films in a multilayered system, sandwiched between two layers of a material with higher PMA and exchange. In this interpretation, each simulation layer can describe a heavy metal/ferromagnet/heavy metal trilayer, which is interconnected to neighboring layers by a combination of direct exchange and stray fields.⁵⁰ Alternatively, but experimentally more challenging, our model could also use a sandwiched bulk polar magnet model with C_nv symmetry, allowing iDMI and hence the stabilization of Néel skyrmions.^14,15 Nevertheless, fabrication of complex multilayer systems requires sophisticated optimization of the stack. Here, we focus on the fundamental properties of a coupled bobber system as a proof of concept to disentangle different magnetic textures.

For simulations in a synthetic antiferromagnetic environment, the 20 lateral cells were divided into two layers of 10 cells. The two layers directly adjacent to the interface are coupled by the Ruderman–Kittel–Kasuya–Yosida (RKKY) interaction, described by the following contribution to the total micromagnetic energy as readily available⁵¹ in the MuMax3 framework⁵² used here:

with σ being the RKKY interaction coefficient, d_NM the thickness of the non-magnetic spacer, which is set to 1 nm in all of our simulations, and $mi⃗$ the normalized magnetization of the top and bottom coupled layers. The RKKY coefficient is set as σ = −2 × 10⁻⁴ J m⁻², lying in the center of a range of values studied for coupled 2D SAF skyrmions.⁴⁰ Note that the method used to simulate the SAF disregards the thickness of the spacer layer in the stray field calculation and is geometry dependent due to interactions of the spin textures with edge modes (see the supplementary material). For a direct comparison to experimental results, it can be included at the expense of smaller cells and longer simulation times, and should model the experimentally investigated sample design.⁵³ To assess the validity of our approach, H–L and L–L bi-bobbers were generated in systems with a gap of 1 nm between the RKKY coupled layers and otherwise the same parameters as stated above. These showed nearly identical bi-bobber energetics as the simulations without gaps.

The magnetization was initialized as Néel-type SAF skyrmion tubes in the chiral layers and as uniform magnetization in the top and bottom high PMA layers, causing the tubes to pinch off and form bobbers. The uniform layers were therefore disregarded for simulations of skyrmion tubes to obtain tubes with a constant diameter throughout the material. Relaxation was performed by integrating an adapted LLG equation, containing only the damping term, using the Runge–Kutta method until a local energy minimum was reached, and a subsequent minimization by the steepest descent method.^52,54

As briefly described before, this work will mainly focus on an antiferromagnetically stacked system of Néel bobbers, and compare their dynamic breathing modes to those of antiferromagnetically coupled Néel skyrmion tubes. Both structures are illustrated in Figs. 1(a) and 1(b), respectively.

The polarity p of skyrmionic excitations is defined as the spin direction of the core. The chirality c is defined by the sense of rotation of the domain wall separating the core from the ferromagnetic background. In SAF environments, the two coupled objects are of opposite polarity, but different combinations of chiralities are possible. As an example, a depiction of p = −1 bobbers with right-handed (c = 1) and left-handed (c = −1) chirality is given in Figs. 1(c) and 1(d), respectively. Experimentally, the sign of c can be chosen in heterostructures by interface engineering.^55,56 An often overlooked effect of skyrmions is their asymmetric stray field, which depends on their chirality. This is illustrated by simulations of the stray field as depicted in Figs. 1(c) and 1(d). The sense of rotation of the latter allows for the stray field to form a closed loop on the top surface, hence lowering the protruding stray field, whereas the former does not. It should be noted that this is reversed at the bottom surface. This effect is well-known on a macroscopic scale, where chiral arrangements of permanent magnets create two opposite surfaces with different stray field strengths. An illustration of this is given in Fig. 1(e). On the nanoscale, such effects have been reported and used for the distinction between Bloch and Néel skyrmions by magnetic force microscopy^47,48 but have received little attention when it comes to the interaction of skyrmions living in different magnetic layers. By inspecting Figs. 1(c) and 1(d), four different combinations of stackings with differing stray field interactions become visible. Two symmetric cases, where either both SAF bobbers are of the type shown in Fig. 1(c) and have high invasive stray fields at the interface (H–H), or both are of the type in Fig. 1(d) with low stray fields (L–L). Also, two asymmetric cases with mixed fields are possible, whereas those are just mirror images of one another, i.e., only one needs to be studied (e.g., H–L). This distinction of coupled systems also holds for the simulated skyrmion tubes.

Cross-sections of the simulated layouts are shown in Figs. 1(a) and 1(b). A detailed description of material parameters is given in the modeling section. The top and bottom uniform layers in the bi-bobber sample design have higher PMA and exchange stiffness. They are ferromagnetically coupled to their adjacent chiral layers, while the two chiral blocks are antiferromagnetically coupled through a scaled⁵¹ RKKY-type interaction at their interface [see Eq. (1)]. The ultrathin non-magnetic spacer defining the interface is not simulated. For simulations of the coupled skyrmion tubes, the uniform layers at the very top and bottom were not simulated, as shown in Fig. 1(b).

To assess the robustness of SAF skyrmions and bi-bobbers, the three different stacking types were first compared to each other in terms of their stability in external fields. This presents the easiest experimental means to judge their properties. Interestingly, the field required to destroy the objects does highly depend on the chirality, or rather the combination of chiralities. Both, SAF bi-bobbers and skyrmion tubes, were the least stable when stacked in the H–H configuration. For the specific set of material parameters chosen for the bi-bobber, even very low external fields of 1 mT are enough to transform them into a single bobber in one layer and a uniform state in the other, i.e., delete one-half of the bi-bobber. For such low fields, the distinction as to which layer transforms to the ferromagnetic state is random. Simulations without stray fields result in stable excitations of the same size in the H–H configuration, confirming that the previously described relaxation to a single bobber is indeed governed by the stray field. Also, in the symmetric cases H–H and L–L, the spins on the domain walls of the skyrmionic objects in the top and bottom layers are not antiferromagnetically aligned. This leads to a higher contribution to the RKKY energy compared to the asymmetric H–L case. However, only in the H–H configuration, strong stray fields of both bobbers (skyrmion tubes) are directly exposed to each other, resembling two repelling magnetic poles, indicating that the symmetric interfacial stray fields destabilize the spin texture. This can be seen even more clearly for the H–H combination of skyrmion tubes, where energy minimization leads to a lateral shift of the coupled skyrmions with respect to one another. In an attempt to obtain non-displaced asymmetric coupled skyrmion tubes, simulations were conducted with higher σ = 1 ⋅ 10⁻³ J m⁻² and different skyrmion tube sizes to modify the relative strength of RKKY-coupling and dipole–dipole interaction. However, the stray field always spatially separated the skyrmion tubes. Therefore, this configuration is omitted in the following analysis.

On the other hand, H–L and L–L combinations were stable for a range of different external fields. A comparison is given in Fig. 2. The bobber diameters in each SAF layer were evaluated at the interface where they reach their highest value. In the H–L configuration, SAF bi-bobbers were stable at fields exceeding 500 mT, whereas in the L–L configuration one of the bobbers was destroyed at 380 mT. Generally, SAF bi-bobber diameters were smaller in the L–L configuration compared to the counterpart values at same external fields in the H–L configuration.

Additionally, the asymmetric stray field has another effect that becomes visible in Fig. 2: In the H–L configuration, the bobber sizes at zero field do not coincide. The above also applies for SAF skyrmion tubes.

Furthermore, to compare the stability between SAF skyrmion tubes and bi-bobbers, zero-field simulations of both were carried out for a reasonable range of material parameters. The parameters of the uniformly magnetized top and bottom layers necessary to stabilize SAF bi-bobbers were not altered. Skyrmionic excitations that exceeded a diameter of 70 nm were considered unstable, since the relaxation of these states in bigger disks leads to constantly increasing skyrmion diameters. Stabilization of these bubbles was mainly attributed to the disk edge and that was not accounted for in the stability region plot.

In Fig. 2(b), the stability comparison is shown for the H–L configurations, which were the most stable for both objects, for ranges of A_ex = (1 − 20) pJ m⁻¹, D = (1.5 − 3.5) mJ m⁻² and K_eff = (1 ⋅ 10¹ − 1.1 ⋅ 10³) kJ m⁻³. As can be seen, SAF bi-bobbers approximately have an equal region of stability compared to that of SAF skyrmion tubes within the tested parameter ranges. Due to the proximity of the bi-bobbers to a high exchange and high PMA layer, the stability region of the former shifted toward lower values for A_ex and K_eff for the chiral layers.

As pointed out earlier, a promising approach to distinguish different types of 3D magnetic textures may be via the excitation of breathing modes. In this section, this breathing mode response to a varying out-of-plane external magnetic field B_ac(t) is presented. To quickly and reliably study excitations up to f_max = 100 GHz, a sinc pulse-shape is chosen for the AC magnetic field of the following form: $Bac(t)=Bdc+Bac0⋅sinc[2πfmax⋅(t−t0)]$⁠, where B_dc is a constant external field offset. Before applying a varying field component, the spin textures were allowed to relax in a field of B_dc by the MuMax3 steepest descent minimization method. The field offset guides the initial object size before application of a pulse. Higher offset values will increase the diameter disparity of skyrmionic objects in the antiferromagnetically coupled layers, as already depicted in Fig. 2. Additionally, to eliminate any frequency contribution that might arise due to a relaxation process, a temporal offset of t₀ = 1 ns was added to the sinc pulse. An amplitude of $Bac0=0.5$ mT was used for all simulations.

To measure the frequency response, first, the out-of-plane spin component mz(t) was averaged over the full simulated geometry at every simulation time-step. Also, to distinguish between responses of each layer, the top and bottom chiral layers were averaged separately at every simulation time-step. To resolve modes up to 100 GHz, the chosen time-step fulfills the Nyquist–Shannon theorem τ_step < 1/(2 ⋅ 100 GHz). Next, the averaged changes of the out-of-plane spin components ⟨δm_z(t)⟩ were decomposed into frequency space by Fourier transformation, and the power spectral density (PSD) was computed.

First, the PSD of SAF skyrmion tubes and SAF bi-bobbers in an H–L configuration at fixed external fields are presented for each antiferromagnetically coupled layer in Figs. 3(a)–3(d). As apparent, the resonance modes of SAF bi-bobbers are much richer compared to SAF skyrmion tubes. Since the internal modes of skyrmions depend on their size,⁴⁵ and bobbers contain skyrmions of multiple diameters, this result is in line with our expectations. Figures 3(a) and 3(c) are results at 50 and 20 mT for SAF bi-bobbers and SAF skyrmion tubes, respectively. These are the fields necessary to compensate for the asymmetric stray fields due to the Halbach effect, thus canceling out the size offset of the bobbers or tubes. Hence, with both bobbers and tubes having approximately the same size, their resonances also align. Figures 3(b) and 3(d) show the PSDs at 130 mT. In this configuration, the external field will decrease the magnetic structure size for both bobbers and tubes in the bottom layer, while increasing it in the top layer. As reported for 2D skyrmions, smaller sizes can increase frequencies in a specific band of breathing modes lying around 25 GHz.⁴⁵ In our simulations, the peaks start splitting accordingly. It should be noted that the increasing trend of frequencies for smaller skyrmions is not true for all internal modes. As an example of the peak splitting, a zoomed-in plot of the splitting of peaks 2 and 3 is presented for each layer in the inset of Fig. 3(b). Here, it also becomes clearer that frequency responses in one layer also lead to a small peak in the other. For example, the single bobber in the top layer, as plotted in the inset, would not have an additional peak at around 25 GHz. This can be confirmed when viewing the internal mode in Fig. 3(a). Being antiferromagnetically coupled to a smaller bobber resonating at that value triggers this additional response peak. This effect is mutual and appears in both layers. It also appears for coupled skyrmion tubes as can be very clearly seen in Fig. 3(d) at 8 and 19 GHz.

Next, the internal modes are visualized by a cell-by-cell calculation of the PSD. Color-coded cross sections of the results help to get a better understanding of the resonances in real space. Figures 3(e) and 3(f) show the two modes that appeared for the SAF bi-bobber. The figures only show internal modes in the top layer. A mirrored twin of each internal mode can exist in the bottom layer. Every peak in the PSD can be represented by a combination of these modes. Both modes have strong contributions of the BPs. In the mode shown in Fig. 3(e), the layers at the interface with the biggest skyrmion diameter resonate with a comparable strength to the BP. The contribution decreases across the bobber with diminishing skyrmion sizes. In the mode shown in Fig. 3(f), this behavior is reversed. Additionally, a ring outside the bobber participates in the breathing motion, although to a much lower extent. When viewed in a cross-section, this ring is visible as the two lighter spots in Fig. 3(f). As discussed above, none of the internal modes appeared in a single layer. Every resonance triggered a smaller response in its antiferromagnetically coupled neighbor. To illustrate this effect, the full resonance images are shown for two frequencies in Figs. 3(g) and 3(h). None of the resonances was a mixed mode between types shown in Figs. 3(e) and 3(f). In other words, the enforced resonance was always of the same type. All modes were, to a varying degree, hybridized with edge modes,^24,40,45 indicating that the geometry of the sample will need to be taken into account when using this method to distinguish different magnetic textures.

For a direct comparison, Figs. 3(i)–3(k) show the same visualization for the breathing modes of skyrmion tubes at 130 mT. The frequencies at which these were evaluated are highlighted accordingly in Fig. 3(d).

Finally, a full set of PSDs for external fields ranging from 200 mT to 200 mT is depicted as 2D color-coded maps in Fig. 4. With this, dispersion relations are given for all frequencies up to 60 GHz. Again, all the effects described earlier can be verified. For the SAF bi-bobber frequency dispersion, as shown in Fig. 4(a), modes 3 and 4 cross at around 50 mT. This is the field at which the asymmetric stray field is compensated, and both bobbers have roughly the same size. The triggered frequency response can be clearly seen in the individual PSDs of each layer, presented in Fig. 4(a) on the left. Modes 1b and 1t fully compensate in the combined plot for fields higher than 60 mT. This means, these specific modes can be driven in a strongly coupled breathing phase by external fields. In this range, both layers breathe with the same amplitude and phase. Note that the symmetry of the cell-by-cell PSD shows the primary components of the oscillations to be breathing modes; however, small contributions produced by gyration might still be present, though undetectable.

When comparing the dispersion of SAF bi-bobber and SAF skyrmion tubes, the strongest bi-bobber resonance modes (3 and 4 in Fig. 4) evolve linearly with the external field, as indicated by the dashed lines for the top and bottom layer in Fig. 4(a) on the left-hand side. Although the skyrmion tube has only two modes in the studied range, the dispersion is more complex. There appears to be no linear relationship. Furthermore, at −90 and 110 mT, the system is driven into an unstable state. In the usual case of similarly sized tubes, the RKKY energy is minimized by both tubes remaining aligned on top of each other. Here, with the size discrepancy between the two tubes becoming considerably larger. The smaller skyrmion tube can now move within the area of the bigger tube without significant changes in RKKY energy. However, the stray field energy is minimized for the smaller skyrmion tube being as far away from the center of its bigger counterpart as possible. Since nothing is breaking the symmetry in the simulations, very small fields and numerical error can be enough to decide in which direction the small tube will shift. Such lateral shifts have been described before for the H–H configuration. In this state, low frequency gyration modes appear, which are not studied further in this work and may require additional extended investigation. It seems that the frequencies of modes 1 and 2 in this configuration gradually decrease toward zero at −90 and 110 mT, respectively, right before entering the unstable state. Although at higher external fields this strong disparity between the sizes and displacement is still holding, stable breathing modes reappear.

Additional to previous considerations, Fig. 4 contains simulations for the L–L configuration shown for the SAF bi-bobber and SAF skyrmion tubes in Figs. 4(b) and 4(d), respectively. Due to symmetric stray fields, bobbers or tubes in both SAF layers have the same size at zero-field, and the previously described crossing points shift to 0 mT. Generally, the number of modes is the same for both configurations. Comparing the SAF bi-bobber responses, the L–L configuration modes are shifted toward lower frequencies. Apart from modes 3 and 4, also modes 1 and 2 have a linear response and a crossing point. Other than in the H–L configuration, there is no full compensation of any modes in the studied parameter ranges. SAF skyrmion tube responses do not show the same instability in the L–L case. Frequency minima are reached at lower values of external fields (−50 and 50 mT) and do not decrease toward unstable states. As previously described, this effect is related to the stray field, which is lower in this configuration.

In this work, we investigated a coupled bi-bobber state in a synthetic antiferromagnetic material that was sandwiched between two uniform ferromagnetic layers. The resulting bobber configurations exhibit strongly chirality-dependent stray fields that cause the coupling to become dependent on the winding orientation of the underlying bobbers, thus affecting the spin texture’s stability. This offers a unique platform for the study of collective excitations. We further report that breathing mode excitations provide a valid approach for distinguishing between 3D magnetic textures that otherwise appear very similar.

See the supplementary material for the effects of edge modes and information on the dispersion of saturated states.

P.B. acknowledges financial support from EPSRC (UK) through a DTG. Funding by the International Max Planck Research School for Science and Technology of Nano-Systems and the Max Planck Center for Quantum Materials is gratefully acknowledged. This project was partially funded by the Deutsche Forschungsgemeinschaft [DFG, German Research Foundation; Project No. 403505322, Priority Programme (SPP) 2137]. The authors thank Professor Dr. Jonathan Leliaert (Ghent University) for fruitful discussions and feedback on the model and code, and for proofreading the article.

The authors have no conflicts to disclose.

P. Bassirian: Formal analysis (equal); Investigation (equal); Methodology (equal); Writing – original draft (equal); Writing – review & editing (equal). T. Hesjedal: Supervision (equal); Writing – original draft (supporting). S. S. P. Parkin: Supervision (equal); Writing – original draft (supporting). K. Litzius: Project administration (equal); Supervision (equal); Writing – original draft (equal); Writing – review & editing (equal).

The source scripts for bi-bobber breathing can be found publicly available at https//doi.org/10.5281/zenodo.6956156 to be used as reference. The full set of analyzed data that support the findings of this study are available from the corresponding author upon reasonable request.