Magneto-active elastic shells with tunable buckling strength

Shell buckling is central in many biological structures and advanced functional materials, even if, traditionally, this elastic instability has been regarded as a catastrophic phenomenon to be avoided for engineering structures. Either way, predicting critical buckling conditions remains a long-standing challenge. The subcritical nature of shell buckling imparts extreme sensitivity to material and geometric imperfections. Consequently, measured critical loads are inevitably lower than classic theoretical predictions. Here, we present a robust mechanism to dynamically tune the buckling strength of shells, exploiting the coupling between mechanics and magnetism. Our experiments on pressurized spherical shells made of a hard-magnetic elastomer demonstrate the tunability of their buckling pressure via magnetic actuation. We develop a theoretical model for thin magnetic elastic shells, which rationalizes the underlying mechanism, in excellent agreement with experiments. A dimensionless magneto-elastic buckling number is recognized as the key governing parameter, combining the geometric, mechanical, and magnetic properties of the system.

Shells are curved thin structures that can withstand extreme loading conditions due to the interplay between bending and stretching deformation 1,2 through the so-called 'shell effect' 3 . Thin shells are an ubiquitous structural element in engineering 4 and also widely observed in nature across length scales, from viruses 5 , and capsules 6,7,8 , to pollen grains 9 , and plants 10 . While curvature is the key ingredient underlying the excellent mechanical performance of shells, it is also responsible for the catastrophic (subcritical) nature of their elastic instabilities 11 . Consequently, shells are high sensitivity to imperfections 1,12,13,14 . For over a century, the pressure buckling of a spherical shell has been a long-standing canonical problem of elastic (in)stability 1,2 . When the in-out pressure differential across a shell exceeds a threshold value, the shell loses its load-carrying capacity 1 . The ensuing collapse is unpredictable, occurring at loads that are significantly lower than classic theoretical predictions 14 . Consequently, measuring and predicting the critical buckling pressure has proven to be nontrivial due to the high imperfection sensitivity.
In 1915, Zoelly 15 derived the theoretical prediction for the critical buckling pressure of a perfect spherical shell, of radius R and thickness h, loaded under a uniform pressure p: where E and ν are the Young's modulus and Poisson's ratio of the material, respectively. Notwithstanding this classic result, given the high imperfection sensitivity, buckling pressures measured in experiments, p max , have long been found 1,12,13,14 to be much lower than the prediction from Eq. (1). This discrepancy generated a long debate in the shell mechanics community that lasted for nearly four decades until the disagreement between theory and experiments was finally attributed to imperfections 16 . Given this intrinsic mismatch, it became standard to define an empirical quantity, the so-called knockdown factor, which is always smaller than unity, spreading a wide range 14,17 : 0.05 ≤ κ d ≤ 0.9. * To whom correspondence should be sent: pedro.reis@epfl.ch In engineering, significant efforts have focused on improving predictions for the knockdown factor and understanding how it is affected by imperfections 1,12,13,14 . A breakthrough came only recently, from experiments, when the combination of a rapid prototyping technique for spherical shells 18 and its adaptation to seed precisely designed defects led to a quantitative relationship between the knockdown factor and defect geometry 17 . This study demonstrated that if imperfections can be measured precisely, then the knockdown factor can be precisely predicted using an appropriate shell theory 19,20,21,22,23,24,25 , thus opening the door for less conservative designs of shell structures. In parallel, nondestructive probing methods have recently been proposed to experimentally access the stability landscape of shells 26,27,28 , while accepting the inevitable presence of multiple imperfections to set the knockdown factor. Furthermore, bilayer shells undergoing differential swelling were recently shown to have a varying knockdown factor during a transient, demonstrating how the knockdown factor can be modified post-fabrication, albeit not reversibly on-demand 29 . Still, to date, the knockdown factor is regarded as an intrinsic structural property dictated by imperfections imparted during fabrication or encoded at the design stage.
Here, we propose a robust mechanism to dynamically tune the buckling strength of shells (i.e., their knockdown factor), by coupling elastic deformation and magnetic actuation, to provide active control and liberate the predestination of imperfections. Our framework leverages recent advances in magneto-elastic mechanical systems -including shape-programmable materials 30,31,32,33 , biomedical devices 34 , and soft robots 35,36,36,37 -, which we augment in the context of shell mechanics. Specifically, we manufacture spherical shells made of a hard-magnetic elastomer (hereon referred to as magnetic shells) and characterize their critical buckling pressure under an applied magnetic field. We demonstrate that the knockdown factor can be modified externally by adjusting the magnitude and polarity of the field. To rationalize our results, we develop a magnetic shell theory, providing a physical interpretation of the interplay between shell mechanics and magnetic fields. Furthermore, we uncover the magneto-elastic buckling number (a dimensionless quantity that combines the magnetic, elastic, and geometrical properties), which acts as the single governing parameter of the system.

Experiments on tuning the knockdown factor of pressurized magnetic shells
In our experiments, we position a hemispherical shell of radius, R = 25.4 mm, and thickness, h = 321.2 µm, in between a set of Helmholtz coils (Fig. 1a,b). These two coils impose a steady, uniaxial and uniform magnetic field (flux density vector B a = B a e 3 ) on the shell, perpendicularly to its equatorial plane (see Methods and Supplementary Notes S1 and S2 for details). The shell is made of a magnetorheological elastomer (MRE), a composite of hard-magnetic NdPrFeB particles (volume fraction c = 7%) and vinylpolysiloxane (VPS) polymer (see Methods). The full 3D configuration and the corresponding cross-section profile of the shell are visualized through X-ray micro-computed tomography (µCT, 100 Scanco Medical AG), a representative example of which is presented in Fig. 1c. To measure the buckling strength, we depressurize the shell using a pneumatic-loading system under imposed-volume conditions and measure the associated pressure sustained by the shell (Methods and Supplementary Note S2). Figure 1d presents the load-carrying behavior of the shell, characterized by the pressure (p) as a function of the volume change (∆V ), both of which are normalized, respectively, by the classic buckling prediction, p c , and the corresponding volume change immediately prior to buckling 19,20 , ∆V c = 2π(1 − ν)R 2 h/ 3(1 − ν 2 ). The onset of buckling corresponds to the maximum of each curve, p max = p max /p c , and the accompanying pressure drop indicates the loss of load-carrying capacity of the shell.
In the absence of magnetic field (B a = 0 mT), the representative shell shown in Fig. 1c buckles at the dimensionless pressure level of p max = 0.44. This value is far below the classical prediction of 1 (p max = p c ), the reason being that we intentionally seed a precisely engineered dimple-like defect at the pole during manufacturing (Methods and Supplementary Note S1), so as to consider the influence of imperfections in a controllable manner. Hence, the geometry of the shell deviates slightly from a perfect hemisphere by an amplitude of δ over the region of polar angle 0 ≤ ϕ ≤ ϕ • . The half angular width of the defect is then measured by ϕ • , usually rescaled as ϕ • = ϕ • ( 12(1 − ν 2 )R/h) 1 2 17,19,21,25 . For this test shell, the defect geometry is characterized as δ = δ/h = 0.39 and ϕ • = 11.7 • (ϕ • = 3.2) using an optical profilometer (Methods and Supplementary Note S1). Although the introduced defect is too small to be seen with the naked eye, the buckling strength of the shell is dramatically reduced, by more than a half due to the high sensitivity to imperfections 1,17,19 .
We proceed by focusing on the effect of the magnetic field on the buckling instability. The magnetized shell possesses a residual flux density of B r = B r e 3 (B r = 63 mT, see Methods and Supplementary Note S1), and is responsive to an external magnetic field, with a significant modification of the buckling pressure (Fig. 1d) compared to the no-field case. When the field vector is parallel to the axis of magnetization (i.e., B a and B r are in the same direction), we observe an increase of the critical loads by 12% and 24% under the flux densities B a = 33 mT and B a = 66 mT, respectively. Meanwhile, the accompanying pressure drop at the onset of buckling is gradually reduced, eventually disappearing for B a = 66 mT. Therefore, the shell is strengthened by the applied field, and the buckling event becomes less catastrophic. By contrast, for the opposite field polarity (i.e., B a and d (1) (2)  (2) a set of Helmholtz coils and depressurized using a pneumatic-loading system, comprising a (3) syringe pump and (4) a pressure sensor . The Helmholtz coils are driven by the current I in the direction shown in b, output from a DC power supply (5). b, Computed field of the vector (arrow) and magnitude (contour line) of the magnetic flux densityB a generated by the coils, under current I = 1 A (Supplementary Note S2). The field is uniaxial and uniform in the central region with a flux density of B a = B a e 3 . Given axisymmetry, the field is represented in the x-z plane (y = 0). Coordinates x and z are normalized by the coil radius R c = 46.5 mm. c, Photograph of the representative imperfect shell specimen (top left). Image of its cross-section (top right) scanned through x-ray µCT and the reconstructed 3D view (bottom, one quarter was artificially hidden to aid with visualization). Scale bars are 5 mm. d, Loading curves of pressure, p = p/p c , versus volume change, ∆V = ∆V /∆V c , for the shell shown in c. The peak value represented by the open circle on each curve is the critical buckling pressure. The buckling test is performed at different levels of external flux density B a . Inset: Schematic of a shell subjected to a magnetic field, which is opposite (left) or parallel (right) to the shell magnetization.
B r are in opposite directions), the shell is weakened; the critical load decreases with a consequent more abrupt pressure drop past the buckling event (B a = {−33, −66} mT in Fig. 1d). These findings demonstrate that the intrinsic buckling strength of a magnetically active shell can be modified (increased or decreased), under an external magnetic field, on-demand.
To further explore the effect of magnetism on the buckling strength of pressurized spherical shells, we test shells with different defect geometries, over a range of external flux densities. As shown in the photographs of the specimens in Fig. 2a, the defect amplitudes are varied as δ = {0.21, 0.27, 0.39, 1.26, 2.61} during manufacturing (Methods and Supplementary Note S1). The axisymmetric defect profiles, w, are presented in Fig. 2b, defined as the radial deviation of the measured shell profile from a perfect hemisphere. Figure 2c presents the corresponding knockdown factor measurements, κ d = p max , as a function of the external flux density, B a . Naturally, due to imperfection sensitivity, the shells exhibit distinct knockdown factors for different values of δ. However, in the presence of the magnetic field, we consistently observe an increasing or decreasing knockdown factor over the explored range of δ. Within the range of B a accessible in our experiments, κ d can be changed up to approximately ±30%, with respect to the no-field case. These experimental results demonstrate that the knockdown factor of a shell can be dynamically tuned, as an extrinsic quantity, by adjusting the polarity and strength of the applied magnetic field, via a robust mechanism that is insensitive to imperfections.
profile of the shell with δ = 1.26 and ϕ • = 3.2. c, Knockdown factor, κ d , versus flux density of the applied field, B a , for shells containing the defects shown in b. Symbols and lines represent experimental results and theoretical predictions from the axisymmetric shell model, respectively. The error bars in the experimental data correspond to the accuracy of pressure measurements and the standard deviation of the measurements of shell thickness and Young's modulus; the error bands of theoretical predictions correspond to the standard deviation of the measurements of defect amplitude and shell thickness. Insets: Schematics of directions of the shell rotation vector, q, and the magnetic torque vector, τ = −kq, near the pole, for the cases of B a opposite (left) or parallel (right) to B r .

Theory of axisymmetric hard-magnetic shells
To rationalize the experimental results presented above, we develop a theoretical model that predicts the response of hard-magnetic axisymmetric thin shells under a combination of mechanical and magnetic loading. We consider the Helmholtz free energy of ideal hard-magnetic soft materials 35,38 , comprising an elastic energy term (related to material deformation) and a magnetic energy term (describing the work to align the residual magnetic vector along the external magnetic field). The shell model is developed by reducing this threedimensional energy to the 1D profile curve of the middle surface of the shell (detailed derivation provided in Supplementary Note S3). The dimensional reduction of the Kirchhoff-Saint Venant elastic energy 39 , valid for small strains and large displacements, with the potential of live pressure was reported recently 40 . In the present study, we focus on the magnetic energy term, which, for 3D scale-free materials, is written as 35,38,41 where F is the deformation gradient 39 . We describe the axisymmetric shell profile by the coordinates (ρ,z),ρ being the radial coordinate in the plane (x, y) perpendicular to the axis of axisymmetry (z). Accented quantities (e.g.,å) refer to the undeformed configuration of the shell. The reduced 1D magnetic energy, normalized by πEhR 2 /(4(1 − ν 2 )), can be derived as (Supplementary Note S3) whereå is the area measure, and F is a dimensionless function that depends on the initial and deformed configurations of the shell. From Eq. (3), we identify a magneto-elastic parameter which represents the intrinsic magneto-elastic coupling of the system. Equilibrium equations can be generated by minimizing the total energy for all the possible displacements of the shell, which were solved via the Newton-Raphson method (see Methods and Supplementary Note S3.1.4) 40 .
It is important to highlight that, in our simulations of the 1D model presented above, we use the geometric and physical parameters of the system measured in experiments (see Methods), without any fitting parameters. The defect profile in the region 0 ≤ ϕ ≤ ϕ • is described analytically by w/h = −δ 1 − ϕ 2 /ϕ 2 • 2 ( Fig. 2b), derived based on a simple plate model (Supplementary Note S1). Excellent agreement is found between theory and experiments (Fig. 2c). The variation of the knockdown factor for shells with different defect geometries under the magnetic field is accurately predicted by our shell model, which is, therefore, able to describe the intricate coupling between elasticity, magnetism, and the nonlinear mechanics of thin shells.
Physical interpretation of the reduced magnetic energy Even though our theoretical model can predict the buckling strength of shells, the reduced magnetic energy is highly nonlinear, and the mechanism underlying the change of knockdown factors still needs to be clarified. Therefore, we proceed to expand the integrand of the reduced magnetic energy density in Eq. (3), F(ρ, ϕ ,z, ϕ , ρ, ϕ , z, ϕ ), up to second order in the displacement field (Supplementary Note S3). By examining the role of each term in the buckling instability, we conclude that only the following second-order term dominates the buckling of the shell, where q is the rotation vector of a material fiber of the shell 2,42 , and τ = −k(ϕ)q can be interpreted as a distributed (dimensionless) torque applied by the external magnetic field. This torque is a linear function of the shell rotation with a deformation-independent pre-factor k = λ mχ (ϕ), which we can interpret as the (dimensionless) stiffness of distributed rotational springs, whereχ(ϕ) =ρ, 2 ϕ /(ρ, 2 ϕ +z, 2 ϕ ).
Under pressure loading, material fibers tend to rotate, thereby increasing the magnetic torque (Supplementary Note S3). Whether the torque reacts to restore the undeformed orientation of the material fibers or to rotate them further away from the initial orientation, depends on the sign of the equivalent stiffness k (insects in Fig. 2c). When B a and B r are in the same direction, the stiffness k is positive, thus ensuring that τ is opposite to q, thereby counteracting buckling. As a result, we observe the strengthening of shells with increasing critical loads (cases with B a > 0 in Fig. 2c). Conversely, when B a and B r and in opposite directions, the negative stiffness (k < 0) leads to a torque τ that acts to increase the rotation q. Indeed, in this regime, buckling occurs at lower pressure levels (cases with B a < 0 in Fig. 2c).

Scaling analysis of the change of knockdown factor
With a physical interpretation of the magnetic energy for shells at hand, we now employ scaling arguments to more clearly rationalize how the knockdown factor is modified by the magnetic field for different radius-tothickness ratios. Since the magnetic energy interacts with the live pressure potential to alter the knockdown factor, we balance the two energy terms. The dimensionless magnetic energy can be shown to scale as U m ∼ λ m h/R, while the change of knockdown factor with respect to the non-magnetic case, ∆κ d = κ d − κ d | B a =0 , scales as ∆κ d ∼ ∆p/E(R/h) 2 (Supplementary Note S3). By balancing the scalings for the potential of ∆p, that is U ∆p = ∆p/E, and the magnetic energy U m , we find ∆p/E ∼ λ m h/R, translating into From Eq. (6), we define the magneto-elastic buckling number, which governs the knockdown factor under combined pressure and magnetic loading of magnetic shells. The scaling ∆κ d ∼ Λ m provides a scale-invariant description of how the magnetic field modifies the buckling pressure for shells with different radius-to-thickness ratios.

Robustness of the mechanism to geometric imperfections
We set out to investigate the role of imperfections in the buckling of magnetic shells, which we now address with a more systematic parametric exploration. We fabricate shells over a wide range of the defect amplitude (0.1 ≤ δ ≤ 4) and measure the corresponding knockdown factors, κ d , from buckling tests. In parallel, we run 1D simulations for the same material and geometrical properties. Figure 3a illustrates the relationship between κ d and δ at different levels of external flux density, included in the magneto-elastic buckling number, Λ m , of Eq. (7). The signature of imperfection sensitivity is still observed in the presence of the magnetic field: κ d decreases dramatically when the defect amplitude increases in the regime of relatively small defects (0 < δ < 1). Surprisingly, the results of the change in knockdown factor under the applied magnetic field (∆κ d ) presented in Fig. 3b are significantly less sensitive to defects; ∆κ d becomes approximately constant for δ > 1. This finding suggests that the magnetic interaction between the shell and the external field is nearly unaffected by the intrinsic imperfection sensitivity, ensuring a robust tuning of the knockdown factor.
Moreover, the results in Fig. 3a,b include data for two sets of shells with a different combination of material and geometric properties (E = 1.15 MPa, R/h = 79.1 for MRE-22 shells and E = 1.69 MPa, R/h = 91.3 for MRE-32 shells; see Methods). Still, the κ d and ∆κ d data collapse for these two sets, given that the value of Λ m is the same for both. This collapse supports Λ m as the governing parameter, combining the mechanical, magnetic, and geometrical properties of the system. Throughout the analysis, the theory is in excellent agreement with experiments.
A robust way to quantify the effect of the magnetic field on the change in knockdown factor is to focus on the plateau in Fig. 3b, ∆κ p , defined as the average of ∆κ d over the extent of the plateau-like region. Figure 4a  different radius-to-thickness ratios and defect widths. While the radius-to-thickness ratio strongly influences the change in knockdown factor, the defect width has little effect on it (Supplementary Note 4). In Fig. 4b, we plot the raw data of Fig. 4a as a function of the magneto-elastic buckling number Λ m , finding that the different results for shells with different geometries fall onto a master curve. This collapse demonstrates that the plateau change of knockdown factor is governed by the single parameter Λ m , independently of the geometry of the defect. The superposition of experimental results on the theoretical predictions further validates this collapse. Moreover, consistently with the scaling in Eq. (6), a linear relationship is found between ∆κ p and Λ m , with a slope of 1.26 ± 0.01 obtained via linear fitting. This master curve serves as a concrete design guideline for magnetic shells, summarizing the effect of the magnetic field in tuning the buckling strength of pressurized spherical shells with different material and geometrical properties.
In closing, we have shown that the buckling strength of shells can be dynamically tuned by exploiting the interplay between mechanics and magnetism. The proposed mechanism represents a robust way to gain control on a property of thin shells that has long been regarded as intrinsic, the knockdown factor. By performing precision experiments on hard-magnetic elastomeric shells and developing a theoretical model, we unveiled the essential feature at the base of the mechanism, a distributed torque induced by the magnetic field throughout the shell. Moreover, we showed that a dimensionless quantity, the magneto-elastic buckling number, emerges as the key governing parameter, summarizing the geometrical, mechanical, and magnetic properties of the system. We envision that our mechanism can also be used to gain control on the structural life of a shell, by tuning the magnetic field as critical conditions approach. More generally, we believe that the approach of coupling mechanical deformation and magnetic fields, already successfully applied for beams 32,36,38,43 , films 33 , and soft materials 35 , and employed here to tune the knockdown factor of shells, can be extended to modify the instabilities of other structures such as rods, plates, and non-axisymmetric shells 30,31 , for which research efforts are currently ongoing. We anticipate that magnetic structures with predictable and tunable buckling behavior will enrich the design of advanced materials with new functionalities.

Methods
Preparation of the MRE material. Our experimental samples were fabricated using a MRE material composed of a mixture of hard-magnetic NdPrFeB particles (average size of 5 µm, MQFP- 15-7-20065-089, Magnequench) and Vinylpolysiloxane (VPS, Elite Double, Zhermack), a silicone-based polymer. MREs made of VPS Elite Double 22 or VPS Elite Double 32 are referenced throughout the text as MRE-22 or MRE-32, respectively. We prepared the MRE with the following steps. First, the initially non-magnetized NdPrFeB particles were mixed with the liquid VPS base (1:1 mass ratio) using a centrifugal mixer (ARE-250, Thinky Corporation); for 40 s at 2000 rpm (clockwise), and another 20 s at 2200 rpm (counterclockwise). Secondly, the solution was degassed in a vacuum chamber (absolute pressure below 8 mbar), and then cooled down to room temperature (22 ± 0.4 • C) to avoid any changes of viscosity due to thermal disturbances from the previous steps. Thirdly, VPS catalyst was added to the mixture, with a ratio of 1:1 in weight to the VPS base. After another mixing step for 20 s at 2000 rpm (clockwise) followed by 10 s at 2200 rpm (counterclockwise), the final mixture was ready for sample fabrication. The final mass concentrations of NdPrFeB particles and VPS polymer were 33.3 wt% and 66.7 wt%, respectively. Pouring of the liquid MRE mixture onto the mold during shell fabrication (see below) was done after a set waiting time (100 s for MRE-22 and 20 s for MRE-32), so as to increase the viscosity to the desired value. After the above preparation steps, the curing of the polymer mixture occurred in approximately 20 min at room temperature.
Physical properties of the MRE. The densities of MRE-22 and MRE-32 were measured using a pycnometer to be ρ MRE−22 = 1.61 ± 0.05 g/cm 3 and ρ MRE−32 = 1.63 ± 0.03 g/cm 3 , respectively. The density of the NdPrFeB particles is 7.61 g/cm 3 (provided by the supplier). Under the concentrations of the prepared mixture, the volume fraction of NdPrFeB particles was calculated to be 7.05% for MRE-22 and 7.14% for MRE-32. When saturated, the effective residual flux density of our MRE, B r , was assumed to be the volume-average of the total residual flux density of individual NdPrFeB particles (0.90 T, as reported by the supplier). Given the chosen volume fraction, B r = 63.2 mT for MRE-22 and B r = 64.0 mT for MRE-32. The Young's moduli of MRE-22 and MRE-32 were measured to be E = 1.15 ± 0.04 MPa, and E = 1.69 ± 0.06 MPa, respectively, using a combination of cantilever and tensile tests.
Fabrication of the magneto-active shell specimens. Our shells were fabricated by coating the underside of a flexible negative spherical mold (radius 25.4 mm) with liquid MRE; see the schematic in Supplementary  Fig. S1 and Supplementary Note S1 for details. Whereas our technique was inspired by previous work 17,18 , our molds also contained a soft spot (thin circular region) at the pole to produce a precisely-engineered defect.
With the negative spherical surface of the mold facing down, the gravity-driven viscous flow of the polymer yielded a thin layer of MRE on the mold. Before the MRE fully cured (≈ 13 min after pouring), the mold was depressurized from within using a syringe pump. As a result, the soft spot of the mold deformed to produce an axisymmetric geometric imperfection at the pole of the shell. This defect becomes 'frozen' in the shell upon curing. The amplitude of the defect can be systematically varied by adjusting the pressure imposed on the mold during curing ( Supplementary Fig. S2). The width of the defect can be varied by changing the pillar used in mold manufacture with a different radius (Supplementary Note S1). The amplitude and width of the defect can be set independently. To make the shell magnetically active, we saturated the NdPrFeB particles, already in the solidified MRE, using an applied uniaxial magnetic field (4.4 T, perpendicular to the equatorial plane) generated by an impulse magnetizer (IM-K-010020-A, Magnet-Physik Dr. Steingroever GmbH).
Geometry of the shells. The radius of the shells (R = 25.4 mm) is set by the mold used during fabrication. Their thickness was characterized using a microscope (VHX-5000, Keyence Corporation) after cutting off narrow strips (≈ 2 × 6 mm 2 ) near the pole. The average measured thickness was h = 321.2 ± 5.1 µm and h = 278.2 ± 2.8 µm for the MRE-22 and the MRE-32 shells, respectively. Hence, the corresponding radius-to-thickness ratios were R/h = 79.1 and R/h = 91.3. The 3D profile of the outer surface of each fabricated imperfect shell was characterized using an optical profilometer (VR-3200, Keyence Corporation). The geometry of the defect was computed as the radial distance between the measured profile and the corresponding best-fit spherical surface. Specifically, the amplitude, δ, and half angular width, ϕ • , of the defect were determined by fitting the analytical description, , to the experimental measurements (additional details provided in Supplementary Note S1). The ranges of geometric parameters for our specimens are as follows: 0.14 ≤ δ ≤ 3.2 and ϕ • = 11.7 ± 0.1 • (ϕ • = 3.2) for the MRE-22 shells; and 0.17 ≤ δ ≤ 3.4 and ϕ • = 11.6 ± 0.1 • (ϕ • = 3.4) for the MRE-32 shells.
Generation of the magnetic field. A uniaxial uniform magnetic field was generated in the central region of two identical customized multi-turn circular coils (square cross-section, inner diameter 72 mm, outer diameter 114 mm, and height 21 mm), configured in the Helmholtz configuration. The coils were set concentrically, with a centre-to-centre axial distance of 46.5 mm, and connected in series, powered by a DC power supply providing a maximum current/power of 25 A/1.5 kW (EA-PSI 9200-25 T, EA-Elektro-Automatik GmbH). The flux density of the magnetic field was varied in the range −66 mT ≤ B a ≤ 66 mT by adjusting the current output of the power supply, from −10 A to 10 A. The flux density was measured using a Teslameter (FH 55, Magnet-Physik Dr. Steingroever GmbH), along both the axial and radial directions. In parallel to the experiments, we simulated the generated field using the Magnetic Fields interface embedded in the commercial package COMSOL Multiphysics (v. 5.2, COMSOL Inc.), based on Ampère's Law (see Supplementary Note S2 for details). Excellent agreement was found between experiments and simulations ( Supplementary Fig. S3).
Buckling tests. To quantify the critical buckling conditions of our magnetic shells, we positioned each shell in between the Helmholtz coils under a steady magnetic field, and then depressurized it under prescribed volume conditions by a pneumatic loading system (see Supplementary Note S2 for details). The applied pressure was monitored by a pressure sensor (785-HSCDRRN002NDAA5, Honeywell International Inc.). The level of depressurization was increased, under a steady magnetic field, until the shell buckled. The critical buckling pressure was defined as the maximum value of the loading curve (pressure vs. imposed volume change).
Theory and numerics. The reduced 1D energy U, representing the sum of the elastic and magnetic energies of the system, together with the potential of the live pressure, is obtained by means of dimensional reduction (see Supplementary Note S3 for details). Equilibrium equations are generated by minimizing the total energy for all possible displacements of the shell, and solved via the Newton-Raphson method using the commercial package COMSOL Multiphysics (v. 5.

S1 Supplementary Notes on the Experimental Protocols
In this section, we elaborate on the protocols that we have developed to fabricate and characterize our hardmagnetic shells. In Section S1.1, we present the coating technique used to fabricate the shell specimens. Then, in Section S1.2, we discuss the methodology followed to characterize the geometry of the engineered defects of the fabricated shells.

S1.1 Fabrication protocol
We fabricated magnetic shell specimens using a customized coating protocol (see schematic in Fig. S1) that was inspired from our previous work S1, S2 , while also including some important new modifications. The protocol contains 3 steps -(i) Fabrication of the mold, (ii) Fabrication of the shell, and (iii) Magnetization of the shell -that are detailed next.
(i) Fabrication of the mold: An elastic mold was manufactured as the negative of a rigid hemisphere (radius 25.4 mm, stainless steel ball, TIS-GmbH) by pouring liquid VPS solution (base to catalyst ratio of 1:1 in weight, fully mixed and degassed) into a boxed-shaped container with the hemisphere placed on its bottom side (Fig. S1a). A cylindrical pillar (radius 4.41 mm) was 3D-printed and positioned at a distance of ≈ 0.5 mm above the north pole of the steel ball, in order to produce a circular thin region in the mold. After curing of the VPS, the mold was detached from both the ball and the pillar (Fig. S1b). Then, we activated the surface of the mold using air plasma generated by a plasma cleaner (PDC-002-HPCE, Harrick Plasma) for 2 minutes and treated it with TMCS (Trimethylchlorosilane, Sigma-Aldrich) in a vacuum chamber for 1 hour. This treatment reduced the adhesion strength between the mold and the MRE to facilitate the demolding of shell specimens during step (ii).
(ii) Fabrication of the shell: The prepared mixture of MRE was poured into the negative spherical mold produced during step (i); see Fig. S1c and the Methods section of the main article. The pouring was done after waiting for a set time, counting from the completion of mixture preparation (100 s for MRE-22 and 20 s for MRE-32). The purpose of this waiting time was to increase the viscosity to a desired value S2 . With the spherical surface facing down, the gravity-driven viscous flow yielded a thin layer of MRE on the mold (Fig. S1d). Before the MRE was fully cured, at about 13 min after pouring, we applied constant negative pressure to the mold by extracting the air inside using a syringe pump (Fig. S1e). Under this depressurization, the soft spot (thin region) of the mold deformed, thereby producing an axisymmetric geometric imperfection in the shell at the pole. This defect was 'frozen' in the shell after curing. During this depressurization stage, the bulk regions of the mold (other than the soft spot) exhibited negligible deformation. Shells made of MRE-22 or MRE-32 were fabricated using molds made out of the corresponding VPS polymer.
(iii) Magnetization of the shell: After fabrication, the shells were then magnetized permanently by saturating the NdPrFeB particles in the MRE using an applied uniaxial magnetic field (4.4 T, perpendicular to the equatorial plane, Fig. S1f ) generated by an impulse magnetizer (IM-K-010020-A, Magnet-Physik Dr. Steingroever GmbH). To constrain the shell in its natural undeformed configuration, the process was conducted before demolding by placing the cured shell (still adhered to the pressurized mold) inside the magnetization coil (MF-AsA-φ70×120mm, Magnet-Physik Dr. Steingroever GmbH). We then released the pressure loading and injected VPS mixture into the shell from its bottom to make a thick band (≈ 6 mm in thickness, Fig. S1g). This latter step ensured clamped boundary condition at the equator during buckling tests. Finally, the magnetic hemispherical shell containing a geometric defect was peeled off from the mold (Fig. S1h).
Fig. S1: Schematic diagram of the fabrication protocol for precisely imperfect magnetic shells. a, Fabrication of a negative spherical mold containing a soft spot (thin region) at its pole. b, Unmolding from the metal ball. c, The mold is filled with the liquid mixture of MRE. d, The mold is turned upside-down, which leads to a gravity-driven flow of the MRE that generates a thin coating on its surface. e, The soft spot (circular, thin, plate-like region) of the mold deforms under the applied negative pressure, p δ , to produce the geometric imperfection in the shell. f, The cured shell, while still adhered to the pressurized mold, is magnetized permanently using an applied magnetic field. g, A thick base is made at the equator by injection of additional VPS polymer from the base of the mold. h, The final magnetic shell containing a geometric defect is peeled off from the mold. The relevant geometric quantities are shown.
Systematic variation of the defect geometry. The fabrication protocol presented above allows us to independently vary both the defect width (by changing the radius of the pillar used in the mold manufacture; Fig. S1a) and the defect amplitude (by adjusting the pressure imposed on the mold; Fig. S1e). For the fabrication of the MRE-22 shells, we applied pressure in the range of 1 ≤ p δ ≤ 9.75 kPa, which yielded defect amplitudes ranging from δ = 0.14 to δ = 3.2 (Fig. S2).  Throughout the study, the half angular width of the defects was fixed at ϕ • = 11.7 ± 0.1 • by setting the radius of the soft spot of the mold at 4.41 mm (corresponding to the target half angular width of 10.0 • set by the radius of the pillar). The minor difference between the obtained value for half angular width and the targeted value can be attributed to small deformations during fabrication near the boundary of the spot, which was not ideally clamped by the surrounding soft bulk regions of the mold. For the MRE-32 shells, we varied the defect amplitude in the range 0.17 ≤ δ ≤ 3.4 by changing the pressure in the range 2 ≤ p δ ≤ 13 kPa, with a fixed half angular width of ϕ • = 11.6 ± 0.1 • . An empirical linear relationship was found between the defect amplitude and the applied pressure (the slope depends on the bending stiffness of the thin region of the mold). These results serve as a convenient set of design guidelines to fabricate the specimens.

S1.2 Characterization of the defect geometry
Measurement of the defect profile. To characterize the defect geometry of each shell specimen, the 3D profile of the shell (outer surface) was measured using an optical profilometer (VR-3200, Keyence Corporation). A perfect spherical surface was fitted to the measured data in the region away from the defect located at the pole of the shell. The defect profile, w δ , was defined as the radial distance between the fitting spherical surface and the measured profile. Due to the axisymmetry of the shell, the 3D profile was averaged latitude-wise to obtain the 2D profile.
A simplified plate model to describe the generated defect geometry. We used a simple model to provide an analytical description of the geometry of the defect. Note that, during pressurization of the mold, the bulk region outside of the soft spot exhibits negligible deformation compared to the soft spot. Moreover, the spherical cap of this soft spot is sufficiently small (with respect to the radius of curvature of the shell), such that it can be regarded as a (nearly) flat plate. Hence, we isolate this thin region and model it as a boundary-clamped circular flat plate loaded by a uniform pressure. The defect profile is determined by the deflection of this plate under the pressure loading. The governing equation can be obtained from classic Kirchhoff-Love theory for an isotropic and homogeneous plate under pure bending S3 , where w δ is the out-of-plane displacement of the plate, and x,y are the two in-plane coordinates, respectively. Solving Eq. (S1) for the specific case of a circular plate of radius, r • , and thickness, t, under uniform pressure loading, p δ , and clamped boundary conditions, yields the following expression for the plate deflection S3 where D = Et 3 /12 1 − ν 2 is the bending modulus of the plate, and r = x 2 + y 2 is the radial coordinate. Note that we are using t to designate the thickness of the circular plate, so as not to confuse it with the thickness of shells, h, used elsewhere in the study. Recognizing the maximum deflection at the center of the circular plate, w δ (0) = −p δ r 4 • /(64D) as the defect amplitude δ, and r/r • as ≈ ϕ/ϕ • (from geometry), the profile of the dimple-like defect, for ϕ ≤ ϕ • , can be then written as where ϕ • is the half angular width of the defect, and ϕ is the angular distance (polar angle) measured from the center of the defect. Outside of the defect (ϕ > ϕ • ), we have w δ = 0. We determined the amplitude, δ, and half angular width, ϕ • , of the defect by fitting Eq. (S3) to the experimentally measured defect profiles, finding excellent agreement between the two (see Fig. 2b of the main article). From Eq. (S2), we also note that δ = w δ (0) scales linearly with p δ , which allows us to vary the defect amplitude δ during shell fabrication by adjusting the pressure applied to the mold. In Fig. S2, we present the values of δ measured after demolding versus p δ for a series of MRE-22 and MRE-32 shells, confirming this linear relationship. The slope depends on the bending stiffness of the thin region of the mold. The non-zero intercept on the horizontal axis is due to the natural deflection of the thin region, which might be induced by the pre-stresses developed in the mold during the curing of VPS polymer.
So far, we have provided a simplified theoretical model to aid in the characterization of the profile of the defects in our imperfect shells. This analysis provides an excellent descriptor for the functional form of the defect geometry, even if the underlying two parameters (amplitude and width of the defect) are determined through fitting. In addition to aiding in setting the design parameters for fabrication, the geometric profile of the imperfect shells provided by this model is used directly as an input for the shell model (Section S3). In Section S4, we will investigate how the geometric profile affects the onset of buckling by independently varying the underlying geometric parameters in the numerics.

S2 Supplementary Notes on the Experiments
In Section S2.1, we provide information on the generation of the magnetic field in our experiments, including the coils setup with the simulations and experiments performed to characterize the generated field. Then, in Section S2.2, we detail the apparatus and procedure for the buckling tests, with a focus on the generation of the magnetic field.

S2.1 Generation and characterization of the magnetic field
Description of the Helmholtz coils used to generate the magnetic field. Two customized multi-turn circular coils with a square cross-section were set up in the Helmholtz configuration to generate a uniaxial uniform magnetic field in their central region. Each coil was manufactured by winding an aluminum circular spool with an insulated magnet wire: 18 turns in the axial direction and 19 turns in the radial direction (Repelec Moteurs S.A.). The magnet wire (enamelled wire, Isomet AG) had a circular cross-section of diameter 1 mm for the copper core and thickness 0.094 mm for the outer insulation layer. The final dimensions of the coil were: 72 mm in inner diameter, 114 mm in outer diameter, and 21 mm in height. The two identical coils were set concentrically, spaced with an axial centre-to-centre distance of 46.5 mm. They were connected in series and powered by a DC power supply that provided a maximum current/power of 25 A/1.5 kW (EA-z 9200-25 T, EA-Elektro-Automatik GmbH). The strength of magnetic field was varied by adjusting the current output from the power supply.
Numerical simulation of the generated magnetic field. We used the commercial package COMSOL Multiphysics (v5.2, COMSOL Inc.) to simulate and characterize the magnetic field generated by the coils described above. These simulations were conducted using the Magnetic Fields interface embedded in COMSOL, based on Ampère's Law. Given the axisymmetry of the system, we modeled the 2D cross-section of the two coils, encompassed by a circular region of air (relative permeability 1) that was significantly larger than the size of the system (radius ten times larger than that of the coils). In this setting, we only consider the copper parts (relative permeability 0.999994) with current flowing through; i.e., the aluminum housing of the wires, with a relative permeability of ≈ 1, were neglected. Since the strength of the magnetic field is linearly proportional to the current, I, we normalized the simulations by setting I = 1 A.
Experimental characterization of the generated magnetic field. Experimental measurements of the field strength were carried out to characterize the field generated in the experiments, as well as to validate the simulations. An axial Hall probe (HS-AGB5-4805, Magnet-Physik Dr. Steingroever GmbH) was moved along the central axis (vertical) of the two coils (inset of Fig. S3a), and the corresponding flux density was measured by a Teslameter (FH 55, Magnet-Physik Dr. Steingroever GmbH). Similarly, a transverse Hall probe (HS-TGB5-104010, Magnet-Physik Dr. Steingroever GmbH) was used to measure the flux density along the radial (horizontal) direction (inset of Fig. S3a). These measurements agree well with the simulation results, along both the axial and radial directions, as shown in Fig. S3.   S3: Comparison between the experimentally measured and the computed flux density of the generated magnetic field. The experimental data (symbols) and simulation results obtained using COMSOL (solid lines) along the axial centre line (x = y = 0) and the radial centre line (y = z = 0) are presented in the plots a and b, respectively. A current of 1 A is applied to the coils. The radius of the coils is R c .

S2.2 Buckling tests under combined magnetic and pressure loading
Experimental apparatus for the buckling tests. We measured the critical buckling load of the fabricated magnetic shells under pressure and magnetic loadings by positioning the shell in between the Helmholtz coils, with its pole located in the central region of the field (see Fig. 1a,b). The pneumatic loading system comprised a syringe pump (NE-1000, New Era Pump Systems Inc.) and a differential pressure sensor (total error band ±3.7 Pa, 785-HSCDRRN002NDAA5, Honeywell International Inc.), which were both connected to the shell.
Protocol of the buckling test. For each experimental run, under a given external magnetic flux density, the syringe pump was set to extract the air inside the system at a constant flow rate of 0.4 mL/min, so as to depressurize the shell until it buckled. During depressurization, the pressure differential between the outside (atmospheric pressure) and inside of the shell was monitored by the pressure sensor at an acquisition rate of 8 Hz. For each tested shell, we obtained the evolution of pressure loading as a function of the imposed volume extracted by the syringe. The critical buckling pressure is defined as the peak value of this signal, representing the load-carrying capacity of the shell.
Avoiding overheating of the coils during the experimental tests was crucial. Otherwise, the resulting increase of temperature in the apparatus would have caused thermal expansion of air in the pneumatic loading system and compromised the results. This issue was particularly important for the stronger shells (at higher values of B a ). Noting that the flow rate of 0.4 mL/min for the pressurization was set constant for all experiments, these experimental could run for up to 70 s, which was far too long for the coils to be switched on. This increased running time for some of these runs meant that applying the strong magnetic field continuously throughout the full duration of the test was not feasible. Fortunately, auxiliary experiments demonstrated that the measured buckling pressure under the presence of the magnetic field was insensitive to when the field was switched on along the loading path, as long as this was done prior to buckling. Hence, to obtain reliable results, the field could be switched off from the start of the experimental test and only switched on in the proximity of p max , prior to reaching this value. To do so, we considered the cases of B a > 0 (stronger shells with respect to no-field) and B a < 0 (weaker shells with respect to no-field) differently. When the shell is strengthened by the magnetic field (tests for external flux density B a > 0 in Figs. 2c and 3 of the main article), buckling occurs at increasing pressure levels with B a . To avoid overheating of the coils, the magnetic field was only applied when the pressure reached the level that was 10 Pa below the critical buckling pressure of the non-magnetic case (B a = 0). On the other hand, for shells weakened by the magnetic field (in the range −66 < B a [mT] < 0 in Figs. 2c and 3 of the main article), we followed a slightly different procedure while still attempting to minimize overheating. Specifically, we first measured the critical pressure for the lowest value considered for the applied magnetic field (B a = −66 mT), corresponding to the weakest conditions (lowest buckling pressure), while switching on the field at the start of the test. For subsequent experiments, we systematically increased the applied field in steps of 6 mT, over the range −66 < B a [mT] < 0 (11 experiments in total). For the i th experiment, we had the field switched off until the pressure reached the value of (p max (i − 1) − 5) [Pa], where p max (i − 1) was the critical pressure measured in the previous (i − 1) th experiment.

S3 Supplementary Notes on the Theoretical Model
This section is dedicated to the theoretical model that we developed to describe the onset of buckling of axisymmetric hard-magnetic shells. First, we describe the kinematics of axisymmetric shells, followed by the dimensional reduction of the underlying Helmholtz free energy (Section S3.1). To gain insight into this nonlinear magnetic energy, we perform an expansion for small displacements and interpret the physical meaning of the relevant terms at the different orders, namely: magnetic pressure, in-plane force, and torque (Section S3.2). A scaling analysis on these terms is presented in Section S3.3 to help rationalize how the knockdown factor changes with the applied magnetic field. In Section S3.4, we examine their role in the load-carrying behavior of magnetic shells. Emphasis is given to the second-order term of the magnetic energy representing the distributed magnetic torque relevant to the shell rotation, which is found to dominate the critical buckling conditions of the shells.

S3.1 Axisymmetric hard-magnetic shell model
Dimensional reduction of the Helmholtz free energy for axisymmetric shells. As a starting point for the derivation of our hard-magnetic shell model, we use the Helmholtz free energy for ideal hard-magnetic soft materials S4,S5 . This specialization of the free energy comprises an elastic part related to material deformation (hereon referred to as elastic energy) and a magnetic part describing the work to align the residual magnetic moment along the external magnetic field (hereon referred to as magnetic energy). Given that the saturated MRE can be assumed to have the permeability of vacuum S6,S5 , the magnetic potential energy in the reference configuration of volume V can be simplified as the work to align the residual magnetic moment of the material with the external magnetic field. For 3D scale-free materials, this magnetic potential energy is S6,S4,S5 where µ 0 is the permeability of vacuum, F is the deformation gradient tensor, B r is the vector of residual magnetic flux density of the material in the reference configuration, and B a is the vector of external magnetic flux density.
In the following, we seek to derive the magnetic energy counterpart for axisymmetric shells by reducing the 3D description of Eq. (S4) to the 1D profile curve of the middle surface of the shell; a procedure that is commonly referred to as dimensional reduction S7,S8,S9 . This dimensional reduction will be performed based on the Kirchhoff-Love assumptions for thin shells S7 (material fibers normal to the middle surface remain normal upon deformation without stretching nor shrinking), within the framework of small strains and large displacements/rotations. We will then follow a variational approach to minimize the total 1D energy via a weak form implementation, similarly to what we have recently proposed for elastic shells subjected to combined pressure and force-indentation loading S10 .

S3.1.1 Kinematics of axisymmetric shells
Defining the spherical coordinates (ϕ, θ, η) that represent the polar angle, azimuthal angle, and radial distance, respectively, we parametrize the undeformed reference configuration of the shell as a solid body of revolution X(ϕ, θ, η) =r(ϕ, θ) + ηn(ϕ, θ) = (ρ(ϕ) cos θ,ρ(ϕ) sin θ,z(ϕ)) + ηn(ϕ, θ) , (S5) wherer(ϕ, θ) = (ρ(ϕ) cos θ,ρ(ϕ) sin θ,z(ϕ)) represents the undeformed middle surface, with (ρ,z) being the 2D Cartesian coordinates of its profile curve, andn(ϕ, θ) is the undeformed unit normal vector. The covariant base vectors corresponding to X can be derived as Under the Kirchhoff-Love kinematic assumptions, the deformed configuration of the shell can be represented by the current middle surface r(ϕ, θ) = (ρ(ϕ) cos θ, ρ(ϕ) sin θ, z(ϕ)), together with its normal vector n(ϕ, θ), as x(ϕ, θ, ρ) = r(ϕ, θ) + ρn(ϕ, θ) = (ρ(ϕ) cos θ, ρ(ϕ) sin θ, z(ϕ)) + ηn(ϕ, θ) . (S7) The corresponding covariant base vectors are Now, we recall that the deformation gradient is a two-point tensor defined as S11 where G i , defined by G i · G j = δ i j (i, j = 1, 2, 3), are the contravariant base vectors corresponding to X, and G ij are the contravariant components of the metric in the undeformed configuration of the shell. To compute G ij , we first compute the covariant components G ij = G i · G j as resulting in a diagonal metric of the undeformed shell configuration. Hence, the contravariant components can be readily computed as the inverse of the nonzero diagonal covariant components as with all the other components being zero. With these expressions at hand, we will expand the deformation gradient along the thickness coordinate so as to perform the dimensional reduction by integration of the magnetic energy along the thickness. Before doing so, we first need to briefly introduce the geometry of the middle surface, which will be later required for the development of our theoretical model.
The corresponding contravariant base vectorså α are defined byå α ·å β = δ α β (α, β = 1, 2). By definition S12 , the covariant components of the first and second fundamental forms of the undeformed middle surface can be expressed aså αβ =å α ·å β andb αβ = −n, α ·å β , respectively. Moreover, we denote the square root of the determinant of the covariant metric as • a = det ( • a αβ ), which will be used in the area measure. Similarly, if we evaluate the parametrization of the deformed three-dimensional shell, x, at η = 0, we obtain the parametrization of the deformed middle surface, that is r(ϕ, θ) = (ρ(ϕ) cos θ, ρ(ϕ) sin θ, z(ϕ)). The covariant base vectors can be expressed as a 1 = ∂r ∂ϕ = (ρ ,ϕ cos θ, ρ ,ϕ sin θ, z ,ϕ ) , The contravariant base vectors, first and second fundamental forms associated to the deformed middle surface can be computed with the same expressions we showed for the undeformed middle surface, employing the quantities without the accents. Having completed the geometric description of our system, we can now proceed to perform the dimensional reduction of the magnetic energy.

S3.1.3 Dimensional reduction of the magnetic energy for axisymmetric shells
We seek to perform a dimensional reduction of the magnetic potential energy for 3D scale-free materials in Eq. (S4), from 3D to 1D. To do so, consider the specific case for the configuration of the magnetic field in our experiments reported in the main article, where First of all, the elementary volume can be written as where we used the expression of the square root of the determinant of the undeformed three-dimensional metricG for a stack of surfaces S7,S13 , namelyG =å(1 − 2ηH + η 2K ), withH andK denoting the mean and the Gaussian curvatures of the undeformed mid-surface, respectively. Then, expressing the deformation gradient by making use of Eq. (S14), (S10) and (S11), we can write the integrand of Eq. (S4), that is, the density of the magnetic potential, as In the next step of the dimensional reduction procedure, we integrate the magnetic energy along the thickness of the reference configuration of the shell, by making use of Eqs. (S15-S16). The zeroth order terms in η yield terms that are akin to a stretching energy (i.e., proportional to the thickness). The linear terms in η vanish upon integration along the thickness direction from η = −h/2 to η = h/2. The higher-order terms, O(η 2 ), would give bending-like contributions to the energy, which we evaluated numerically to be much smaller than the stretching-like contribution, with no effect on the solution. Then, the integrand (Eq. (S16)) becomes Finally, we obtain the 1D reduced magnetic energy for an axisymmetric shell which we have normalized by πEhR 2 /(4(1−ν 2 )). For simplicity, we have introduced the following dimensionless function in Eq. (S18) which depends on the initial and deformed geometries of the shell. Also, note thatå = ρ 2 (ρ, 2 ϕ +z, 2 ϕ ) is the area measure (i.e., the square root of the determinant of the metric of the undeformed middle surface). A dimensionless magneto-elastic parameter, emerges from this analysis, representing the magneto-elastic coupling for thin magnetic shells.

S3.1.4 Elastic energy and live pressure potential
So far, we have reduced the dimension of the magnetic potential from 3D to 1D, leaving the elastic energy and the live pressure potential aside as these are already known in the shell literature S7,S14,S15,S10 . Next, we recall and briefly report these energies, while emphasizing the underlying assumptions at their origin. The elastic energy U e , as proposed by Koiter S14,S7 , is valid for small strains since it stems from the Kirchhoff-Saint Venant three-dimensional strain energy. Normalized by πEhR 2 /(4(1 − ν 2 )), it reads S10 where the first term is the (dimensionless) stretching energy, and the second term is the (dimensionless) bending energy. Note that we have introduced the trace operator in the undeformed metric "tr " such that, for example, tr (a) = • a αβ a αβ .
The live pressure potential, normalized by πEhR 2 /(4(1 − ν 2 )), can be expressed as S10 Finally, the total dimensionless energy is the sum of the elastic energy, the live pressure potential, and the magnetic energy: Numerical implementation: In our computer simulations, the equilibrium equations can be obtained by minimizing the total energy U for all possible displacements of the shell. These equilibrium equations are solved via the Newton-Raphson method using the commercial package COMSOL Multiphysics (v. 5.2, COMSOL Inc.). Details of the implementation of our numerical procedure are provided in Ref. S10 , albeit without the magnetic contribution derived above.

S3.2 Physical interpretation of the reduced magnetic energy
Above, we derived an expression for the reduced magnetic energy, U m in Eq. (S18). We now seek to provide a physical interpretation of the various terms involved in U m , focusing in the limit of small displacements. Due to axisymmetry, we express the displacement field of the shell profile in the covariant basis {å 1 ,n}, introduced before, as with u 2 = 0, because of axisymmetric deformations. Notice that u α are the contravariant components of the in-plane displacement field (along the covariant base vectorså α defined above), whereas w is the normal displacement S7 . Then, since r =r + u, we can expand the integrand of Eq. (S18) up to second order in the displacement field to find where the constant term, which does not depend on the displacement field, will be neglected from hereon. We will now analyze each term in Eq. (S25), separately, with the goal of understanding their physical significance. Specifically, in the subsequent sub-sections, we will arrive at the following findings (i) The linear term proportional toå 1 · u, 1 , defined as represents the combined action of distributed tangential and normal (dead pressure-like) forces; (ii) The second-order term proportional to (n · u, 1 ) 2 , defined as is equivalent to a distribution of linear rotational springs, and dominates over all the other terms; (iii) The second-order term proportional to (n · u, 1 ) (å 1 · u, 1 ), which we define as represents a distribution of torques that are linear functions of the displacement field.  We start by considering the first-order term in Eq. (S25), defined as U (1) m in Eq. (S26), which is proportional toå 1 · u, 1 . In words, this term represents the component of the longitudinal derivative of the displacement vector along the first covariant base vector (again longitude). Note that, using the covariant basis, we can establish the following identitẙ where we used the Weingarten formulae S12 to expressn, 1 = −b 1 1å1 , withb 1 1 denoting the mixed component along the latitudinal coordinate of the second fundamental form of the undeformed middle surface. We already considered the orthogonality of the spherical coordinates and the fact that the coordinate lines are principal. Then, the approximation for small displacements of the first-order term of the magnetic energy can be written as where we highlighted the two different parts of this term proportional to the in-plane displacement and the normal displacement, respectively, and where is a known dimensionless function depending on the undeformed configuration of the shell.
Magnetic pressure. If we compare the term U (1w) m , proportional to the normal displacement, to the potential of a dead pressure, namely U dead p hå dϕ (positive pressure for a shell under compression), we can define a dimensionless contribution to the magnetic loading (normalized by E) that is analogous to a dead pressure, For a perfect spherical shell withρ = R sin ϕ,z = R cos ϕ,b 11 = −R, andå = R 2 sin ϕ, Eq. (S32) reduces to This equivalent magnetic dead pressure varies across the shell, being zero at the north pole. Notice that the magnitude of this normal force is dictated by the magneto-elastic parameter, λ m , defined in Eq. (S20), with a spatial variation imposed viaψ(ϕ).
Distributed, tangential magnetic forces. Regarding the term proportional to the in-plane displacement, that is U (1u) m , we rewrite it as where we have used the Gauss-Weingarten formula S12 with ∇ 1 denoting the covariant derivative alongå 1 , andΓ 1 11 is a Christoffel symbol of the second kind S12 . Then we integrate Eq. (S34) by parts to obtain The first addend on the RHS of Eq. (S36) vanishes because of the clamped boundary condition at the equator (u 1 | ϕ=π/2 = 0). Consequently, Eq. (S36) becomes which we can also write as where we have defined the distributed tangential forces induced by the magnetic field as In the case of a perfect spherical shell, this term is equal to In conclusion, the first-order contribution to the magnetic energy, can be regarded as a combination of distributed normal and tangential forces that vary across the shell. We now move forward to analyze the second-order terms of U m .
The quantity above is the first covariant component of the rotation vector q = (δå α ·n)å α = −(δn ·å α )å α , which is defined, for example in Ref. S7,S16 , such that the rotation of a material fiber defined by the unit vector t is given by q · t. The rotation vector for small displacement can be expressed as S7,S16 whereb β α is the mixed components of the second fundamental form of the undeformed middle surface. In the case of axisymmetry, we have with q 2 = 0. Therefore, we can write where is a dimensionless geometric function depending exclusively on the undeformed configuration of the shell, similarly toψ(ϕ). In conclusion, the term U (2)A m can be interpreted as the energy potential of a distribution of linear torques τ along the shell, defined as where k, defined in Eq. (S46), is the equivalent rotational stiffness. For a perfect hemispherical shell, k is maximum at the pole and decreases monotonically to zero at the equator.
Magnetic torque proportional to the in-plane displacement. Finally, we turn into the second-order term proportional to (n · u, 1 ) (å 1 · u, 1 ), that is U (2)B m , which can be written as where m(u) is a distributed torque, which is a linear function of the displacement field, expressed as Therefore, the term U (2)B m corresponds to a distribution of torques across the shell, which depend linearly on the displacement field. We will analyze the magnitude of the first and second-order terms studied in this section and show how the term U (2)A m , corresponding to a distribution of linear rotational springs, dominates over the other terms.

S3.3 Scaling analysis
Now that we have a full interpretation of the magnetic energy -Eqs. (S41), (S45) and (S49) -, we perform a scaling analysis to understand how the buckling pressure depends on the magnetic field, the geometry of the shell and its material properties. We will proceed by balancing the potential associated to the live pressure and the magnetic energy.
First, we note that the undeformed surface coordinates, as well as their derivatives, scale as the characteristic radius of the shell; e.g., • ρ ∼ R and • ρ, ϕ ∼ R. Furthermore, as we are interested in the onset of buckling that takes place within the linear deformation regime, we assume that both the in-plane and normal displacements scale linearly with the shell thickness, |ǔ| ∼ h and w ∼ h, all the way up to the onset of buckling.
Scaling of the first-order term, U (1) m . Considering the above assumptions for the displacement field, we find the following scaling for the first-order term of the reduced magnetic energy in in Eq. (S41): noticing a linear proportionality to the thickness-to-radius ratio.
Scaling of the second-order term, U A m . We have to first focus on the scaling of the rotation vector q in Eq. (S43). Given that the majority of the deformation of the shell takes place in the neighborhood of the north pole, within a length that scales with the boundary layer ∼ √ Rh, angles will scale with the angular width of that boundary layer as ∼ h/R. This result means that |q| ∼ h/R or, equivalently, q · q ∼ h/R. Finally, we can derive the scaling of U which scales exactly as the first-order term U m .
Scaling of the second-order term, U Scaling of the potential associated to the live pressure, U p . Since the change in volume scales as R 2 h, the dimensionless change in volume will be of order 1. Therefore, the live pressure potential scales as Scaling law of the change of knockdown factor, ∆κ d , with respect to the non-magnetic case. We can now balance the live pressure potential with the reduced magnetic energy to obtain the scaling for the change of knockdown factor due to the applied magnetic field. Combining the scaling results for the first and second-order terms of the magnetic energy, we conclude that Recalling Eq. (1) of the main article, we note that the classic prediction for the buckling pressure of a perfect spherical shells scales as p c ∼ E(h/R) 2 . Hence, we find that the knockdown factor (defined as κ d = p max /p c , where p max actual critical buckling pressure of the imperfect shell) scales as κ d ∼ p/E(R/h) 2 . Consequently, the change of knockdown factor with respect to the non-magnetic case, κ d = κ d − κ d (B a = 0), scales as ∆κ d ∼ ∆p/E(R/h) 2 . Balancing the scalings for the potential of ∆p, U ∆p = ∆p/E, and the magnetic energy in Hence, the scaling of the change of knockdown factor is From the result in Eq. (S57), we define the magneto-elastic buckling parameter, which governs the knockdown factor under combined pressure and magnetic loading of spherical magnetic shells. The scaling ∆κ d ∼ Λ m allows us to rationalize how the magnetic field modifies the buckling pressure for shells with different radius-to-thickness ratios.

S3.4 Quantification of the effect of the magnetic field on shell buckling
Towards gaining a better physical understanding of the magnitude of the first and second order terms of the reduced magnetic energy, we ran two separate sets of simulations for shells containing a defect of amplitude δ/h = 0.39 and angular width ϕ • = 3.2. First, we used (i) the fully nonlinear magnetic energy (U m , Eq. (S18)). Second, we considered its terms at the different orders: (ii) up to second-order (U In Fig. S4, we plot the loading paths of magnetic shells in both the planep − ∆V (Fig. S4a) and in the planep − w (Fig. S4b), for three different values of the applied magnetic field B a = {66, 0, −66} mT (red, black, and blue lines, respectively). In both plots, the black curves correspond to the classic case of pressure buckling of the shell (B a = 0 mT) S2 . The solid curves correspond to the results from the fully nonlinear magnetic energy U m . When the magnetic field is applied, the loading paths are modified with respect to the zero field case; there is a marked difference between the red/blue curves and the black curves. Consequently, there is a finite change of the knockdown factor, as reported in detail in Figs. 3 and 4 of the main article. The dotted lines correspond to the loading paths obtained considering only the linear term of the reduced magnetic energy; clearly showing that the linear term alone is unable to reproduce the observed finite change in knockdown factor. By contrast, we find that the second-order magnetic energy term U (2)A m is able to accurately reproduce the fully nonlinear behavior; the dotted-dashed curves are close to the solid curves, especially in the pre-buckling regime, up to the onset of the instability. Therefore, the second-order term in the magnetic energy, which includes the contribution of the distributed magnetic torque relevant to the shell rotation, dominates the modification of the critical buckling conditions due to the applied magnetic field.
Distribution of the magnetic torque and its dependence on the applied pressure. The magnetic torque, τ = −k(ϕ)q, is non-uniformly distributed along the shell, depending on the rotation vector, q, and the equivalent rotational stiffness, k(ϕ) = λ mχ (ϕ) (Eq. (S46)). Here, we want to understand how the magnetic torque varies as a function of the rotation vector, for positive and negative applied magnetic fields, so as to shed light on the mechanism that underlies the change in knockdown factor. To do so, we selected the representative shell containing a defect of amplitude δ = 0.39 and width ϕ • = 3.2, which was tested in the experiments (Fig. 1d of the main article).
In Fig. S5, we plot the predictions from our model for the magnitudes of the rotation vector, |q| = √å 11 q 1 , and the torque vector, |τ | = √å 11 τ 1 , along the profile of the shell. Apart from the case of zero magnetic field, two representative values of the magnetic field were considered: B a = 66 mT (Fig. S5a) and B a = −66 mT (Fig. S5b). For each case, we present results for the four pressure levels represented by the open circles marked on the loading curves in Fig. S5c. We find that the shell rotation is localized both at the defect (i.e., at the pole of the shell, ϕ = 0) and near the equator. Due to the non-uniformity of the equivalent rotational stiffness, k, the torque is prominent at the pole but almost zero in the rest of the shell. The magnetic torque increases with the rotation during pressurization.
Note that, for the case when the field is parallel to the shell magnetization (B a = 66 mT), the directions of rotation and torque are opposite. By contrast, for the case when the field is anti-parallel (opposite) to the shell magnetization (B a = −66 mT), the rotation and torque are in the same direction. We have integrated the distributed torque on the shell middle surface and, in Fig. S5d, we plot the normalized resultant torque, τ r = 1 Rh √å 11 τ 1å dθdϕ, as a function of pressure. Note that the initial value of τ r at p/p c = 0 has been subtracted. This initial value represents the total torque in the pressure-free state and was verified to have a negligible effect on the critical loads. Interestingly, we find that the torque does not increase linearly with pressure, with a slope that increases and becomes remarkably steep as the pressure approaches the critical value (p max represented by the dash lines in Fig. S5d), where large rotation occurs.
In conclusion, we find that the change of the knockdown factor is dominantly induced by the torque localized near the defect, while the torque is almost zero throughout the rest of the shell, where rotations are negligible. As such, one should be able to extend the proposed approach, which has been demonstrated on hemispherical shells, for full spherical shells and even shallow shell sections, guaranteed by the localized nature of shell buckling. This would be an extension of classical methods in shell mechanics where, given the shallowness of the buckling mode, studies on the buckling of shallow spherical shells S17 have proven to be successful in predicting the buckling of the corresponding full spherical shells S14 .

S4 Parametric Study of the Effect of the Shell Geometry on Buckling
In this final section, we report the results from a parametric study where we have systematically varied the following two parameters: the shell radius-to-thickness ratio, R/h, and the defect width, ϕ • . We set out to investigate the effect of the shell geometry on buckling of magnetic shells. In Section S4.1, we explored R/h over a wider range ( factor and its change under the applied magnetic field. The effect of the defect width on the critical buckling conditions is also investigated in Section S4.2.

S4.1 Independence of shell radius-to-thickness ratio
Using the 1D model for magnetic shells that we presented in Section S3, we obtained numerical solutions for a series of shells with different values of the radius-to-thickness ratio R/h = {50, 100, 200, 400}, and quantified the predictions for the knockdown factor, κ d , as well as its change with respect to the non-magnetic case, ∆κ d , of the shell. The results are shown in Fig. S6 for different values of the magneto-elastic buckling parameter Λ m = {0.18, 0.09, 0, −0.09, −0.18}. We recall from the definition in Eq. (S58) that Λ m = λ m R/h. We considered defects with increasing amplitudes 0.1 ≤ δ ≤ 4, while fixing the defect width at ϕ • = 3.2. We find that, for any given defect geometry, both κ d and ∆κ d are governed by the parameter Λ m , independently of R/h (there is a near-collapse of the κ d (δ) curves for the different R/h values).

S4.2 Effect of width of the defect on the change of knockdown factor
We also investigated the effect of the defect width, ϕ • , on the critical buckling load using our magnetic shell model. We explored the following values ϕ • = {1.0, 2.0, 3.0, 4.0} and computed the prediction for both the knockdown factor, κ d , of the shell and its change with respect to the non-magnetic case, ∆κ d , for different magnetic loading conditions. In Fig. S7a, we present the predicted κ d (δ) curves as a function of defect amplitude, δ, in the absence of the magnetic field (Λ m = 0). We find that κ d decreases for wider defects (i.e., increasing ϕ • ). The associated change of knockdown factor with respect to the non-magnetic case, ∆κ d , is plotted in Fig. S7b. We find that, when the defect amplitude is larger than a given threshold (determined by the differentiation of ∆κ d with respect to δ smaller than 1% and represented by the symbol on each curve of the plot) the change of knockdown factor becomes independent of the defect width (the curves for the various values of ϕ • collapse). For relatively small defects (δ below the threshold value), ∆κ d depends, but only slightly, on ϕ • . The corresponding variation is lower than 8% for the extreme case explored of ϕ • = 4.0 and Λ m = −0.18; less sensitively in compassion with κ d , which decreases by up to 45% for the same parameter set.