Exploring Bistability in the Cycles of the Solar Dynamo through Global Simulations

The calling card of solar magnetism is the sunspot cycle, during which sunspots regularly reverse their polarity sense every 11 yr. However, a number of more complicated time-dependent behaviors have also been identified. In particular, there are temporal modulations associated with active longitudes and hemispheric asymmetry, when sunspots appear at certain solar longitudes or else in one hemisphere preferentially. So far, a direct link between this asymmetric temporal behavior and the underlying solar dynamo has remained elusive. In this work, we present results from global 3D magnetohydrodynamic simulations, which display both behavior reminiscent of the sunspot cycle (regular polarity reversals and equatorward migration of internal magnetic field) and asymmetric, irregular behavior which we interpret as active longitudes and hemispheric asymmetry in the simulations. The simulations are thus bistable, in that the turbulent convection can stably support two distinct flavors of magnetism at different times, in superposition or with smooth transitions from one state to the other. We discuss this new family of dynamo models in the context of the extensive observations of the Sun’s surface magnetic field with the Solar and Heliospheric Observatory and the Solar Dynamics Observatory, as well as earlier observations of sunspot number and synoptic maps. We suggest that the solar dynamo itself may be bistable in nature, exhibiting two types of temporal behavior in the magnetic field.


Introduction
Since the early observations of the sunspot number, the Sun's magnetic field has been known to follow the fairly regular 11 yr sunspot cycle. Sunspots appear at midlatitudes at the beginning of the cycle, then at latitudes a bit higher during the peak of solar activity, and finally at sites progressively closer to the equator as the magnetic activity wanes (e.g., Hathaway 2011). Polarimetric observations from the National Solar Observatory/Kitt Peak Observatory and magnetograms from the Michelson Doppler Imager (MDI) aboard the Solar and Heliospheric Observatory (SOHO) and the Helioseismic and Magnetic Imager (HMI) aboard the Solar Dynamics Observatory (SDO) have enabled detailed and continuous study of how sunspots appear in pairs of opposite-polarity sense from one cycle to the next, making a 22 yr cycle overall (e.g., Hathaway 2015). In addition, a number of other cycles (less obvious and regular than the 22 yr cycle) have been identified, from a slow ∼100 yr modulation of the peak cycle amplitudes called the Gleissberg Cycle (e.g., Gleissberg 1939;Ogurtsov et al. 2002), to a rapid ∼2 yr periodicity in the global magnetic field (e.g., Ulrich & Tran 2013). On timescales in between, active longitudes (longitudes at which sunspots appear more frequently and with greater strength) persist for several decades (e.g., Henney & Harvey 2002). There are also periods (the longest being ∼50 yr) associated with hemispheric asymmetry, or times at which more sunspots conglomerate in either the northern or southern hemispheres (e.g., Ballester et al. 2005). Finally, there are several recorded grand minima, or times throughout history when solar activity is substantially diminished over protracted intervals of several decades (e.g., Hoyt & Schatten 1998).
All of these complicated features of the solar cycle must have their roots in the interior solar dynamo-the process by which the Sun's interior magnetic field regenerates through dynamical interactions between rotation and convection. Sunspots are believed to be the result of rising toroidal loops of the magnetic field. That their sites of emergence migrate equatorward with the sunspot cycle and that their polarity senses flip every 11 yr suggest an interior toroidal reservoir of magnetism which also migrates equatorward and flips polarity every 11 yr. The additional cycles and behavior mentioned above could also be indicative of corresponding cycles in the interior field.
In the past decade, global 3D magnetohydrodynamic (MHD) simulations of the solar convection zone (CZ) have made significant headway in reproducing aspects of the solar cycle. Using the anelastic spherical harmonic (ASH) code, Brown et al. (2010Brown et al. ( , 2011 have shown it is possible to build strong interior "wreaths" of magnetism amidst chaotic, turbulent flow. These wreaths can be space-filling, with nearly 360°connectivity throughout the spherical shell, and can exhibit regular polarity reversals, during which the longitudinal direction of magnetic field in the wreath flips every few years of simulation time. Also using ASH, Augustson et al. (2015) have achieved wreath-building dynamos which exhibit both equatorward propagation of wreaths at low latitudes and poleward propagation of wreaths at high latitudes-a signature feature of the SOHO/SDO observations (e.g., Hathaway 2015).
Perhaps one of the most significant issues facing global MHD simulations of solar convection is that the wide range of spatial and temporal scales relevant for the fluid makes a direct numerical simulation impossible. Researchers have attempted to address this problem by using various prescriptions for subgrid-scale (SGS) turbulent effects or by restricting the global domain. Warnecke et al. (2014) used the PENCIL CODE to solve the MHD equations in a spherical-wedge geometry, showing that global convection could yield αΩ-type dynamos with the direction of propagation of the interior magnetic field being set by the Parker-Yoshimura rule. Passos & Charbonneau (2014) used the EULAG-MHD code (which incorporates implicit dissipation on the smallest spatial scales to maintain numerical stability) to achieve regular, cyclic polarity reversals in their global millennium simulation. The cycles persisted over very long timescales, with statistical features showing long-term trends reminiscent of the observed solar Gleissberg modulation. Hotta et al. (2016) have explored a large dynamical range for solar convection using their reduced speed of sound technique. They explored the coupling of a near-surface layer of small-scale (<10 Mm) convection and deep, large-scale flows, finding that coherent magnetic structures could persist even in the presence of very small diffusivities.
Another important outstanding issue for solar MHD simulations is that most lack the eruption of the interior magnetic flux and subsequent decay of active regions believed to play an essential role in the global dynamo. In the work of Nelson et al. (2011Nelson et al. ( , 2014, an extension of ASH incorporating a dynamic Smagorinsky treatment of the SGS fluid motions achieved a dynamo in which small loops of magnetism detached from the interior wreaths and rose to the outer surface via magnetic buoyancy instabilities. The polarity, twist, and tilt of the loops displayed statistical properties reminiscent of Joy's law for sunspot emergence. Using a finite-difference spherical anelastic MHD code, Fan & Fang (2014) found that the convection gave rise to super-equipartition magnetic flux bundles that had characteristics similar to emerging active regions on the Sun. More recently, dynamos with wreaths that give rise to buoyant loops have been achieved with Rayleigh (the code used in this work) in modeling the intense magnetism exhibited by M-dwarf stars (Bice & Toomre 2020). These simulations have thus shown it is possible to connect MHD simulations of the interior solar dynamo to the emergence of magnetic flux that is actually observed at the photosphere.
It has generally been difficult for simulations to reproduce the equatorward propagation of magnetic field as seen in the solar butterfly diagram (although see Käpylä et al. 2012 andAugustson et al. 2015 for notable exceptions). Furthermore, two other prominent features of solar magnetic activity-active longitudes and hemispheric asymmetry-have not been systematically explored in simulations. In this work, we report on a new class of dynamo simulations whose magnetism exhibits polarity-reversing cycles with equatorward propagation, as well as a quasi-regular, hemispherically asymmetric cycle, the features of which are suggestive of active longitudes and hemispheric asymmetry. Both cycles are present simultaneously, although usually one cycle is more dominant than the other.
In Section 2, we discuss our numerical methods and the parameter space covered by our simulations, as well as the hydrodynamic progenitor common to each magnetic simulation. In Section 3, we discuss the bistable nature of the magnetism achieved in our dynamo cases. We discuss striking flux-transport (F-T)-like behavior in Section 4. We examine each cycling mode individually in Sections 5 and 6. In Sections 3-6, we refer only to our primary lowest magnetic-Prandtl-number case, and we return to the higher magnetic-Prandtl-number cases in Section 7. In Section 8, we explore the dynamical origins of the two cycles in terms of the production of the magnetic field. In Section 9, we discuss our bistable dynamo simulations in the context of solar observations and present our conclusions in Section 10.

Numerical Experiment
We consider time-dependent, 3D simulations of a rotating, convecting shell of fluid modeled after the solar CZ. We drive the system via a fixed internal heating profile that represents radiative diffusion and injects a solar luminosity (distributed in space) into the layer. The background stratification is marginally stable to convection (i.e., adiabatic), although the heating quickly drives the system to a state of instability, wherein convection carries most of the solar luminosity throughout the bulk of the CZ. The energy injected by the heating is ultimately carried out the top boundary via thermal conduction.
We use the open-source, pseudo-spectral MHD code Rayleigh 0.9.1, which scales efficiently on parallel architectures (Featherstone & Hindman 2016a;Matsui et al. 2016;Featherstone 2018). Our modeling is computationally quite challenging, and we often used of order 10,000 processors on the Pleiades supercomputer. Rayleigh solves the equations of MHD in rotating, convecting spherical shells, expanding the fluid and magnetic variables with spherical harmonics in the horizontal directions and with Chebyshev polynomials in the radial direction. The spherical shell has inner radius r i and outer radius r o . In discussing the mathematics, we use the standard spherical coordinates r, θ, and f-the radial distance, colatitude, and azimuthal angle, respectively-and the corresponding unit vectors e r , q ê , and f ê . Rayleigh implements an anelastic approximation (e.g., Gough 1969;Gilman & Glatzmaier 1981;Jones et al. 2011), which effectively filters out sound waves, making the maximum allowable time step limited by the flow velocity and not the sound speed. The thermodynamic variables are linearized about a temporally steady and spherically symmetric reference state with adiabatic stratification, similar to the stratification of the solar CZ revealed by helioseismology. The pressure, density, temperature, and entropy are denoted by P, ρ, T, and S, respectively. Overbars indicate the fixed reference state, and the absence of overbars indicates deviations from the reference state. Although the reference state is fixed in time, the simulations naturally develop spherically symmetric contributions to the deviations P, ρ, T, and S, which effectively allow the system to achieve a final equilibrium state that is slightly different from the reference state.

Fluid Equations and Boundary Conditions
Under the assumptions of anelastic approximation, the nonlinear fluid MHD equations are given by (e.g., Brown et al. 2010) ) is the magnetic field (also in the rotating frame), c p is the specific heat at constant pressure, and g r ( ) is the gravitational acceleration due to a solar mass M  at the origin. The momentum equation (3) employs the Lantz-Braginsky-Roberts approximation (Lantz 1992;Braginsky & Roberts 1995), which is exact for an adiabatic reference state. We choose the diffusivities n r ( ), k r ( ), and h r ( ) (kinematic viscosity, thermometric conductivity, and electrical resistivity, respectively) to be temporally steady and spherically symmetric, varying with radius as rr 1 2 ( ) . Because it is not computationally possible to resolve the full range of motion down to the Kolmogorov scale, these diffusivities must be regarded as turbulent "eddy" diffusivities. For simplicity, we choose the turbulent n r ( ), k r ( ), and h r ( ) to appear in the MHD equations with the same form as the molecular diffusivities, but with greatly enhanced values. The tensors D and e ij refer to the Newtonian viscous stress and rate of strain, respectively. We use the internal heat source Q(r) to represent heating due to radiation, which is chosen to have the fixed radial profile The choice of normalization constant α ensures that the solar luminosity L  is forced through the domain (see Featherstone & Hindman 2016a).
The equation of state considers small quasistatic thermodynamic perturbations about a perfect gas: γ=5/3 being the ratio of specific heats for a fully ionized gas (i.e., a gas with three translational degrees of freedom). We adopt stress-free and impenetrable boundary conditions on the velocity: For the magnetic field, we match to a potential field at both boundaries: r r r at and , where 0. 8 i o 2

( )
Other magnetic lower boundary conditions, such as a perfect conductor or a purely radial field, could be considered, but there is no obvious choice from a physical point of view.
Finally, we force the entropy perturbation S to satisfy a fixed-flux condition at both boundaries: The outer thermal boundary condition (10) causes a sharp radial gradient in the entropy-i.e., a thermal boundary layerto develop near the outer surface, such that the solar luminosity injected by the internal heating is ultimately carried out via thermal conduction. We note that the microphysics of this thermal boundary layer stand in contrast to the small-scale radiative cooling at the solar photosphere. However, including 3D, nonlocal radiative transfer in global simulations is currently very difficult computationally and thus beyond the scope of this work.

Parameter Space of the Simulation Suite
We discuss three dynamo simulations in the present study, which all start from the same hydrodynamic progenitor. We define W º´--2.87 10 rad s 6 1  as the sidereal Carrington rotation rate, and all our models rotate at three times this rate to ensure a low Rossby number and thus a solar-like differential rotation, in which the equator rotates faster than the poles (e.g., Brown et al. 2010). The transition from solar to antisolar differential rotation as simulations grow increasingly more turbulent has come to be called the "convective conundrum" (O'Mara et al. 2016). The antisolar regime can be avoided by either rotating faster or lowering the luminosity, and we opt for the former. We refer to the hydrodynamic progenitor as case H3 ("H" for hydrodynamic, "3" forW = W 3 0  ). Some input parameters common to case H3 and the three dynamo cases are shown in Table 1.
The three dynamo simulations were initialized from the wellequilibrated hydrodynamic case H3 after about 6400 P rot , using the full MHD equations. These dynamo cases differ in the values of the magnetic Prandtl number n h = Pr m , which takes on the values 1, 2, and 4. For all cases, the kinematic viscosity profile n r ( ) is fixed, effectively keeping the Rayleigh number constant. We set the value of Pr m by varying the value of the magnetic diffusivity at the top of the domain. We label the corresponding dynamo simulations as cases D3-1, D3-2, and Table 1 Input Parameter Values Common to All Four Cases (H3, D3-1, D3-2, and D3-4) r i5 .00 10 10 cm r o6 .59 10 10 cm c p3 .50 10 8 erg K −1 g −1 γ 5/3 r i 0.181 g cm −3 L 3 .85 10 33 erg s −1 M 1 .99 10 33 g Ω 0´-8.61 10 6 rad - ) for the number of radial, latitudinal, and longitudinal grid points, respectively (see Table 2).
Our dynamo cases expand the parameter space explored by Brown et al. (2010)-who examined the case = Pr 0.5 m while varying the rotation rate-to higher magnetic Prandtl numbers. In all cases, the magnetic field was initialized from a random seed field of amplitude ∼1 G ("random" in terms of the amplitude of each individual spectral component).
Solutions to the nonlinear MHD equations in rotating, convecting spherical shells involve highly time-dependent and intricate flow structures. Figure 1 shows a sample of the global flow structure achieved in the hydrodynamic progenitor H3. As seen in the orthographic and Mollweide projections of the radial velocity v r , the flows at any given time consist of Taylor columns between about   20 latitude (also called "Busse columns" and "banana cells" in the literature) and a broken honeycomb network of upflows and downflows at higher latitudes. The Taylor columns transport angular momentum outward (Brun & Toomre 2002;Busse 2002;Matilsky et al. 2019) and drive a strong differential rotation with an associated meridional circulation, as seen in Figures 1(c) and (d).

Bistable Magnetic Cycles in the Dynamo Cases
The magnetism in our dynamo cases exhibits two striking (and distinct) cycling modes. Usually one is dominant, but they can also coexist simultaneously. The phenomenon of bistability in dynamical systems generally refers to the coexistence of two stable equilibrium states, and we use the term somewhat loosely to refer to the two different types of dynamo cycle. In this and the following sections, we discuss case D3-1, which exhibits bistable cycling in the cleanest manner. We return to the higher Pr m cases in Section 7.
Each cycle is associated with a unique morphology of the magnetic field, as seen in Figure 2. In roughly the lower half of the shell at an early time as in Figure 2(a), there is one pair of opposite-polarity wreaths in each hemisphere (a fourfoldwreath structure) confined to within ∼40°latitude of the equator. At this instant, the polarity sense of the wreath structure is largely symmetric across the equator. Each wreath mostly has 360°connectivity, linking the magnetic field in a large torus. The wreaths individually have rms field strengths of ∼5 kG and longitude-averaged rms field strengths ∼2.5 kG.
Later in the simulation, near the middle of the CZ as in Figure 2(b), there is a strong, negative-B f partial wreath (wrapping ∼180°around the sphere), as well as a slightly less prominent positive-B f partial wreath on the opposite side of the sphere. The negative-B f partial wreath is dominant at this time, with a peak strength of ∼40 kG, compared to the positive-B f peak field strength of ∼20 kG. In the longitudinal average, there is a residual á ñ f B that is negative and has a peak magnitude of ∼10 kG. The partial-wreath structure occupies most of the domain in radius but is confined between the equator and 20°south in latitude.
The temporal behavior of the magnetic field as a whole roughly consists of transitions between the two magnetic field Averaged meridional circulation, with the mass flux magnitude in color overlaid on the circulation streamlines. In (d), the red and blue tones correspond to clockwise and counterclockwise circulation, respectively. For any vector field A, we define the meridional part as º + q q A e e A A r r mˆˆ. The angular brackets with the "t" subscript denote a combined temporal and longitudinal average. structures presented in Figure 2. Figure 3(a) shows the evolution in time-latitude space near mid-CZ for á ñ f B (throughout the text, the angular brackets with no subscript denote a longitudinal average at a particular time). After the dynamo has established strong fields from the initial seed field (around = t P 300 rot ), the fourfold-wreath cycle dominates. Each wreath emerges at ∼40°-50°latitude and then migrates steadily equatorward. The newly formed wreaths alternate between positive á ñ f B and negative á ñ f B being dominant, with the time between the appearance of two wreaths of the same polarity-the fourfold-wreath cycle period-equal to ∼25 P rot , or about six months. This is the same time it takes an individual wreath to migrate from its midlatitude starting point to the equator. Effectively, each hemisphere operates on its own (i.e., not necessarily in phase), producing wreaths of a given polarity at midlatitudes that reappear once they reach the equator.
At around = t P 2100 rot , the fourfold-wreath cycles are disturbed by the partial-wreath state, which starts in the south (first arrow in Figure 3(a)). The dominant polarity of the partial-wreath structure reverses quasi-periodically. We can estimate the period (time between two successive states of one dominant polarity) from a visual inspection of Figure 3(a): between t=2950 P rot and t=3350 P rot (an interval slightly larger than the one spanned by the second set of dashed lines), there are roughly five cycles, yielding a period of ∼80 P rot . We note that visual inspections of other intervals would yield slightly different cycle periods, making the period of ∼80 P rot only approximate. Nevertheless, we thus identify a partialwreath cycle with a period about three times as long as that of the fourfold-wreath cycle.
The partial-wreath pair wanders into the north around = t P 2250 rot and flips from north to south once more before the fourfold-wreath cycle becomes dominant for the second time around = t P 4150 rot . The rest of the simulation displays a complex seesaw behavior between the two cycling states. The fourfold-wreath cycle never really disappears, but does significantly decrease in amplitude when the partial-wreath cycle is dominant. Examining Figure 3(a) in detail, it seems that whenever the partial-wreath cycle starts, its field comes originally from the fourfold-wreath cycle. That is, a wreath of one polarity will migrate toward the equator and then significantly grow in amplitude, seeding the following partialwreath cycle with that same polarity dominant. Three instances of this phenomenon are marked by the arrows in Figure 3(a).
Along with the time-latitude panels in Figure 3(a), we show the temporal behavior of the energy in the azimuthal magnetic field in Figure 3(b), averaged over 10% of the shell (by radius) at mid-depth. The fourfold-wreath cycle is seen to correspond to a low-energy state of the magnetic field, with sporadic peaks in energy but no discernible energy cycle associated with the polarity reversals. When the partial-wreath cycle first begins (around = t P 2100 rot ), the field jumps into a high-energy state, with local peaks spaced by roughly half the partial-wreath cycle period.

Polar Caps of Magnetism
Also visible in Figure 3(a) is a weak azimuthal magnetic field at high latitudes on the order of tens of Gauss (see, for example, the positive á ñ f B in the north and negative á ñ f B in the south during the interval from 3500 to 5000 P rot ). These polar caps of magnetism fluctuate in amplitude and occasionally reverse polarity, but not with any regularity and with no discernible relationship to the fourfold-wreath or partial-wreath cycles.
There is a significantly stronger polar cap associated with the meridional magnetic field B m , with magnitude ∼1 kG. Figure 4 shows the radial magnetic field when there are intense polar caps of magnetism. Near the top of the domain (left-hand column), the local structure of B r inherits the honeycomb network of downflows associated with the Taylor columns, while deeper down (right-hand column), the field is more homogeneous. At both depths shown, however, the complicated polar-cap structure on average resembles a dipolar field, with magnetic north lining up with the geometric north. The instantaneous views of B θ (not shown) similarly display a dipolar structure, with B θ positive near the poles in both hemispheres. Figure 5 provides the time-latitude diagram for á ñ B r near the top of the domain during an interval in the middle of the simulation. There is clearly a strong dipolar component to the field at the poles. After a fairly quiescent initial interval, the dipolar field strengthens (at » t P 3750 rot ) a few hundred rotations before the fourfold-wreath cycle becomes dominant  for the second time. The partial wreaths produce flux that is transported to the poles, seen in Figure 5 as intermittent large plumes of red and blue.
Not all partial-wreath cycles produce field that propagates all the way to the poles, and there does not seem to be a one-to-one correspondence between polarity reversals of the partial wreaths and reversals in the polar caps. For nearly every partial-wreath cycle, however, the radial field propagates much farther poleward than the   20 latitude limits that are maintained fairly consistently for á ñ f B . The poleward migration of á ñ B r seen in Figure 5 bears a striking resemblance to the concept the meridional circulation acting as a "conveyor belt" of magnetism in the F-T dynamo framework (e.g., Charbonneau 2014). In that picture, the meridional circulation advects the meridional component of the magnetic field to the poles near the outer surface, where it accumulates and is eventually pumped back down to the equator in the deep CZ, also by the meridional circulation. This poloidal field is then sheared into toroidal field by the differential rotation, starting the cycle anew. In the F-T framework, the meridional circulation time thus sets the cycle period.
However, case D3-1 is more complicated than the simple F-T model for two reasons. First, the meridional circulation time for case D3-1 is on the order of P 400 rot , which is substantially longer than either the partial-wreath or fourfoldwreath cycle periods. Second, although the poles occasionally flip polarity, there is no clean correspondence between the polarity of the field at the poles and the polarity of the field in the equatorial regions during subsequent cycles. Figure 6 shows the evolution of the azimuthal magnetic field B f during one reversal in the fourfold-wreath cycle. At any given time, each wreath has a complicated structure due to stretching and pummeling by the Taylor columns, but generally has the same sense for B f all the way around the sphere, indicating field lines that connect in a large torus. The accompanying snapshots of á ñ f B in the meridional plane indicate that the wreaths are strongest deep in the shell and are confined in radius to the lower half of the CZ and in latitude to within ∼40°of the equator.

Evolution of the Fourfold-wreath Cycle
Each wreath marches steadily equatorward as the cycle progresses. As a wreath of one polarity reaches the equator, a new wreath of that same polarity appears at a higher latitude. In Figure 6(a), for example, there are two opposite-polarity wreaths in the northern hemisphere below ∼30°latitude (roughly equal in amplitude) and a negative-á ñ f B wreath that is just beginning to form above ∼30°latitude. By the time a given wreath migrates from midlatitudes to within ∼10°of the equator, a new wreath of that same polarity has formed at midlatitudes, completing the cycle. Therefore, the configuration in Figure 6(e) has essentially returned to that of Figure 6(a). Figure 7 shows a close-up view (in time-latitude space) of a small time interval containing about eight fourfold-wreath cycles (marked by the first set of dashed lines in Figure 3(a)), including the instants sampled in Figure 6. Here the field is first seen to appear (at its weakest) at roughly 40°latitude north and south and migrates steadily equatorward, attaining its maximum strength between 20°and 30°north and south, while stopping short of the equator at about 10°north and south.
Interestingly, the southern hemisphere is often out of phase with the northern hemisphere, and it is hard to pinpoint whether the fourfold-wreath structure is symmetric about the equator or antisymmetric. Indeed, the simulation exhibits both symmetric and antisymmetric tendencies at different times (for example, compare the symmetric configuration in Figure 2(a) to the antisymmetric configurations shown in Figure 6). At = t P 1200 rot (horizontal center of Figure 7), the wreath structure is clearly largely antisymmetric.
The symmetry of dynamo states in mean-field dynamo theory is usually quantified by the prevalence of the dipole (ℓ=1) and quadrupole (ℓ=2). Our dynamo cases, by contrast, are dominated by higher-order modes that serve to localize each wreath in space, leaving the dipolar and quadrupolar contributions extremely weak. Figure 8 shows the magnetic-field power spectra, averaged over the initial interval during which the fourfold-wreath cycle dominates. The power in longitude-averaged magnetic fields (the m = 0 power) peaks around ℓ=10. The power in the fluctuating magnetic fields ( > m 0 | | ) peaks at the smaller scales, around ℓ=20 for ¢ f B and = ℓ 25 for ¢ B m . The peak around ℓ=10 is consistent with the fourfoldwreath structure as seen in Figures 2(a) and 6: each wreath in a given hemisphere has a latitudinal extent of about 15°, or roughly one-tenth of the 180°in latitude covered by the full sphere. We note that for the longitude-averaged fields, the dipolar and quadrupolar modes ( = ℓ 1 and 2, respectively) are weaker by a factor of ∼10 compared to any of the modes near ℓ=10.
In Figure 9(a), we show the temporal behavior of the power in the longitude-averaged (m = 0) field á ñ f B around ℓ=10 deep in the shell during the fourfold-wreath cycle. Most of the time, the fourfold-wreath structure is symmetric about the equator, with the even modes = ℓ 8, 10, 12 containing up to 60% of the total power for á ñ f B . The two hemispheres are out of phase, however, and the symmetric (odd) and antisymmetric (even) modes can sometimes be roughly equal (e.g., near = t P 800 rot ). Near = t P 1200 rot , the antisymmetric modes dominate, as we noted in discussing Figure 7.
The m=0 power around ℓ=10 is a proxy for the total energy contained in the fourfold-wreath structure. The oscillations in Figure 9(a) thus show that the fourfold-wreath energy oscillates periodically. In Figure 9(b), we decompose the time series for all the power contained around ℓ=10, both symmetric (even) and antisymmetric (odd) modes, into its frequency (or period) components. There is clearly one period that is dominant, which corresponds to an energy cycle of = P P 12.0 energy rot . Because the polarity alternates between energy cycles, the fourfold-wreath cycle period is = = P P P 2 24.0 , 11 fourfold energy rot ( ) similar to the estimate of ∼ P 25 rot mentioned previously.

Evolution of the Partial-wreath Cycle
The partial-wreath component (when present) is strongest at mid-depth. From Figure 3(a), the partial-wreath structure first appears around = t P 2100 rot and cycles between positiveá ñ f B -dominant and negative-á ñ f B -dominant partial-wreath pairs. Unlike the fourfold-wreath structure (which continues cycling regularly even when the partial-wreath structure dominates), the partial-wreath cycling is intermittent, turning off altogether at times, but always returning. There is some variation in the time-latitude behavior of á ñ f B from cycle to cycle, but in general, á ñ f B stays confined to within ∼20°of the equator in a given hemisphere and shows a tendency to propagate toward higher latitudes with time.
Even when the partial-wreath cycling is dominant, the equatorward-propagating fourfold-wreath cycle still persists with its regular period of 24.0 P rot unchanged. However, the fourfold-wreath structure is of much weaker amplitude and is less coherent in the hemisphere where the partial wreaths reside. There is thus a wide dynamical range in the amplitudes of the magnetic field overall. The partial-wreath cycle is associated with longitude-averaged field strengths of ∼10 kG, and when the partial-wreath cycle is dominant, the fourfoldwreath cycle's longitude-averaged field strengths are reduced to several hundred Gauss from ∼2 kG. Figure 10 shows a close-up view of the partial-wreath mode when it is cycling in the south. Here, the propagation of á ñ f B away from the equator is clear. For any given cycle, the partialwreath structure begins roughly at the equator and migrates toward midlatitudes, stopping at ∼  20 . Displaying á ñ f B as in Figures 3(a) and 10 gives a sense of the complex temporal evolution, but it must be complemented by examining the spatial structures of B f before they are averaged. This is particularly true for the partial-wreath cycling, where what appear as polarity reversals in á ñ f B are actually modulations in the relative strengths of the partial wreaths.
To determine in detail how the partial-wreath structure modulates over the course of a cycle, we show in Figure 11 a series of volume renderings of the azimuthal field, in a tracking frame rotating somewhat faster than Ω 0 . In each snapshot, the positive-B f and negative-B f partial wreaths are centered roughly around f=210°and f=30°in longitude, respectively, although "islands" of magnetism associated with each wreath extend around the whole sphere. In the first snapshot (panel a), the positive-B f partial wreath is dominant. Its amplitude and longitudinal extent then steadily decline over the next half-cycle while the negative-B f partial wreath becomes dominant (panels b and c). In the second half of the cycle (panels d and e), the positive-B f wreath expands and strengthens until it is dominant again and the configuration is roughly the same as it was at the beginning of the cycle. In the longitudinal average (e.g., Figure 3(a)), much cancellation occurs, but the sign of á ñ f B still reflects which partial wreath is dominant.
In order to track the partial-wreath structure over multiple cycles, we trace B f with respect to time and longitude in a frame rotating at the same rate as in Figure 11. Figure 12 shows the time-longitude trace for B f at 10°south latitude and at two depths, during an interval when the partial-wreath structure is stronger than the fourfold-wreath structure and cycling in the south. The choice of 10°south latitude cuts the partial wreaths roughly in their core, where the field strengths are highest.
Approximately, the same propagation pattern in timelongitude space occurs both at mid-depth and deep in the shell, despite the field being tracked at the same rate for both depths. This rate is close to the equatorial fluid rotation rate at mid-depth, suggesting that the whole asymmetric wreath structure is "anchored" in the middle of the shell and moves through the fluid as a single entity. At mid-depth, the field pattern is uniformly located at a higher longitude than the field pattern deep down. This displacement does not change with time, suggesting that it is due to the inherent geometry of the wreath rather than the shearing action of the flow.
In the rotating frame chosen, the field structure is approximately stationary. Each partial wreath (either positive B f or negative B f ) is modulated in amplitude over a period of ∼ P 80 rot (the same as the cycle period identified in the timelatitude trace of á ñ f B ), but remains basically at a fixed longitude. Furthermore, the positive-B f and negative-B f wreaths are out of phase with another. At any given time, there are thus two partial wreaths, but in general, when the positive-B f wreath is strong, the negative-B f wreath is weak, and vice versa.
The partial-wreath structure is still not exactly stationary even when tracked in the chosen frame. Figure 12     partial wreaths, implying that the propagation rate of the wreath through the fluid is not constant with time, or that longitudinal extents of each partial wreath shrink and expand throughout the cycle.
In light of Figures 11 and 12, the "polarity reversals" in á ñ f B during the partial-wreath cycle are in fact accomplished by each partial wreath weakening and strengthening mostly in place, with the positive-B f partial wreath oscillating in time  180 out of phase with the negative one. This stands in contrast to the reversals of the fourfold-wreath cycles, which are achieved by full wreaths migrating equatorward from midlatitudes and getting replaced by new wreaths of opposite polarity.

Partial-wreath Cycle Period and Hemispheric Asymmetry
The partial-wreath cycle period is roughly P 80 rot , as identified previously, but it is substantially irregular. We determine the period more formally by calculating the periodicities associated with the relative energy content in the two partial wreaths: º - , 1 2 partial ( ) the "+" and "−" superscripts referring to instantaneous rms averages of the field strength over the regions where B f is positive and negative (respectively), thus measuring the energy content in the positive and negative partial wreaths separately. The mean in the rms is taken over longitude, low latitudes between   25 latitude, and all depths. Unlike the linear longitudinal average á ñ f B , the rms-average A partial is more sensitive to the strongest fields, and thus better determines the relative amplitudes of the partial wreaths as visualized in Figure 11. Figure 13(a) shows the relative energy content over an interval for which no period is easily identifiable in the timelatitude trace of á ñ f B in Figure 3(a), but the partial-wreath cycle is still visible in the well-defined peaks and troughs of A partial . In Figure 13(b), we show the periodogram associated with A partial . There is a distinctive peak in power at º P P 81.5 , 13 partial rot ( ) which we define to be the (dominant) partial-wreath cycle period. Due to the irregularity of the partial-wreath cyclecaused by some cycles being longer or shorter than P partial , as well as significant amplitude and phase modulation-there is a substantial spread in the periods for which the Lomb-Scargle power is large. The dominant period identified in Equation (13) is thus consistent with the previous estimate ∼80 P rot from a visual examination of Figure 3(a). The partial wreaths tend to reside predominantly in one hemisphere or the other, although they are not strictly confined and can in some instances cross the equator. Nonetheless, this phenomenon leads to a substantial hemispheric asymmetry, which we define for our dynamo models as º - the "N" and "S" superscripts referring to instantaneous rms averages of the field strength over the northern and southern hemispheres (respectively). The mean in the rms is taken over longitude, low latitudes (between the equator and 25°latitude for the northern hemisphere and comparably for the southern hemisphere), and all depths. Figure 14(a) shows the temporal behavior of the hemispheric asymmetry for case D3-1 during roughly the latter half of the simulation, when the partial-wreath cycle is usually substantially stronger than the fourfold-wreath cycle. The dominant hemisphere switches between north and south chaotically, with no single period easily seen. We have confirmed this aperiodic nature by computing the periodogram associated with A hem shown in Figure 14(b): there is a wide range of periods that are in the modes ℓ=8, 10, 12 (black curve) and ℓ=7, 9, 11 (red curve). These are the even/odd modes near the peak of the spectrum at » ℓ 10, which contain most of the power. Power is shown during the initial interval for which the fourfold-wreath cycle is dominant, coincident with the top panel in Figure 3.   Figure 3(a). This interval (between 3000 and 3300 P rot ), encompasses about four partial-wreath cycles. The vertical dashed lines denote the instances sampled by the volume renderings of B f in Figure 11. relevant, with no obvious peak in the power spectrum. What is clear from the broad peak near P 2000 rot in the periodogram (and from a visual inspection of Figure 14(a)) is that asymmetry of a single sign-corresponding to partial wreaths cycling preferentially in a single hemisphere-can persist for long timescales, up to ∼1000 P rot .

Bistability Trends with Higher Magnetic Prandtl Number
Because we keep the viscous diffusivity fixed while lowering the magnetic diffusivity, the magnetic structures grow increasingly more complex the higher the magnetic Prandtl number. Figure 15 contains snapshots of the azimuthal magnetic field B f (in Mollweide projection) near the beginning of each of the three dynamo cases. At these early times, the structure of B f consists largely of the fourfold wreaths. As Pr m is increased, the magnetic structures become more shredded and dominated by small-scale features. The local field strength increases, but the global coherence of the field decreases.
Although bistability (in the form of the fourfold and partialwreath cycles) is most clearly evident in case D3-1, similar behavior can be seen in the higher-Pr m cases D3-2 and D3-4. In Figure 16, we show extended time-latitude diagrams for case D3-2, with the time axes and color map scaled to be directly comparable to Figure 3. Case D3-2 exhibits bistability with largely the same characteristics as case D3-1, namely a fourfold-wreath cycle with regular polarity reversals and a uni-hemispherical partial-wreath cycle with irregular polarity reversals. The fourfold-wreath cycle has a weaker signal, however, due to the finer-scale magnetic structures. The system destabilizes after only ∼250 P rot to launch the partial-wreath Figure 11. Volume rendering of B f throughout one partial-wreath cycle, tracked in a frame rotating at 32 nHz faster than the frame rotation rate Ω 0 , approximately equal to the equatorial rotation rate at mid-depth. Only strong fields ( > f B | | 5 kG) are depicted. To show the field everywhere in 3D space and simultaneously emphasize the strong fields, transparency varies linearly with field strength, with structures for which f  B | | 35 kG being completely opaque. The first row (panel a) samples the field when the red partial wreath is dominant, and subsequent samples (panels b-e) are separated by roughly onequarter of the partial-wreath cycle period. The left-hand column shows the full B f profile, and the right-hand column shows only positive B f . The view is from ∼30°north of the equator and 90°longitude. The spherical-shell boundaries in the equatorial plane and the location of  0 longitude are marked by the black curves. cycle, which is less intermittent than in case D3-1, wandering between the northern and southern hemispheres without ever fully shutting down.
From Figure 15(c), the fourfold-wreath structure in case D3-4 is less striking than for the lower magnetic-Prandtl-number cases. Each wreath is also significantly tilted in latitude and thus there is very little signal left after a longitudinal average. The partial wreaths, which appear later in the simulation, are also less coherent. Nonetheless, the fourfold-wreath and partialwreath cycles in case D3-4 can still be detected in sequences of Mollweide projections.
Each of the higher-Pr m cases has substantial polar caps of magnetism, which occasionally reverse polarity but are not clearly tied to either of the two cycles. In time-latitude diagrams for á ñ B r and á ñ q B (not shown), there is a clear tendency for the meridional field to break off from the partial wreaths and move to the poles, just as in Figure 5 for case D3-1.
The biggest global effect of increasing the magnetic Prandtl number is to lower the overall differential rotation while keeping the shape of the rotation profile in the meridional plane largely fixed. Table 2 contains the values of some of the diagnostic nondimensional parameters (as well as the grid resolution and various timescales) characterizing each system. We quantify the strength of the differential rotation as the difference at the outer surface of the rotation rate between the equator and  60 latitude, normalized by the frame rotation rate: . For comparison, for the solar rotation rate determined through helioseismology, DW W = 0.197 0 , if we take W = W 0  and r o to be the radius just below near-surface shear layer (Howe et al. 2000). The magnitude of the differential rotation in H3 is quite strong at 0.228-even greater than the solar value. As Pr m increases, DW W 0 is steadily weakened by the presence of small-scale magnetic structures. These counteract the Reynolds stress produced by the Taylor columns (which tend to drive a strong differential rotation) with feedbacks from the nonaxisymmetric Maxwell stress.
From Table 2, modifying the magnetic Prandtl number has little effect on the resulting level of rotational constraint (parameterized by the Rossby number) or the level of turbulence (parameterized by the Reynolds number). This agrees with the fact that the hydrodynamic structures (in terms of their distribution of length scales and amplitudes) in the dynamo cases are largely similar to those of the hydrodynamic case H3. The main effect of modifying Pr m comes from the degree of complexity in the magnetic structures, as seen by the steady increase of Re m with Pr m .
The bulk energy budget for the simulation suite is shown in Table 3. Case H3 has both the highest kinetic energy from the differential rotation ( f KE ) and the highest kinetic energy overall. For the dynamo cases, as the magnetic Prandtl number is increased, the energy in the meridional circulation (KE m ) and convection (KE c ) stays roughly the same as in case H3, while the energy in the differential rotation ( f KE ) sharply decreases. The magnetic energy density (which is dominated by the fluctuating contribution ME c for all three dynamo cases) steadily increases with increasing Pr m , at the expense of the energy in differential rotation. The energy in the mean magnetic fields f ME and ME m seems to vary in no clear way with Pr m , although the mean energy in the azimuthal field ( f ME ) is significantly less in case D3-4, consistent with the shredded profile of B f seen in Figure 15(c).

The Dynamical Origins of Each Cycling Mode
The real 22 yr sunspot cycle must be the result of a complex interaction between turbulent convection, rotation, and magnetism. Our simulations, although far less turbulent than the Sun, have the advantage that the dynamical origins of each of the two cycling modes identified can be explored in a fair amount of detail. This is because we can probe each term in the dynamical equations with high spatiotemporal resolution and accuracy, something we are unable to do for the Sun. We devote this section to describing the dynamics that lead to the partial-wreath and fourfold-wreath cycles in case D3-1. We focus the discussion around the dominant terms in the induction equation (5) and their physical origins.
We note at the outset that for both cycles, there is very little temporal modulation in either the differential rotation (DR; the main driver of the dynamo through the Ω effect) or the meridional circulation (MC). This stands in contrast to the dynamos of Augustson et al. (2015) and Guerrero et al. (2016), in which significant modulation of the DR (10%-20%) played a significant role. For the fourfold-wreath cycle, the modulation of mean flows is undetectable. For the partial-wreath cycle, there is a ∼5% weaker DR when compared to the fourfoldwreath cycle (due to the stronger magnetic fields), and a ∼1% modulation of the DR from cycle to cycle. The MC develops significant equatorial asymmetry during partial-wreath cycling, but this appears to be mostly a response to the magnetic field and has little effect on the dynamics. Figure 13. (a) Temporal evolution of A partial (Equation (12)) during the interval between 6000 and P 7000 rot in case D3-1. (b) Normalized Lomb-Scargle periodogram of A partial for case D3-1 during the interval from P 5100 rot until the end of the simulation, during which time the partial-wreath pair is consistently stronger than the fourfold-wreath structure.
where "MS," "FS," "TA," and "RD" refer to mean shear, fluctuating shear, total advection, and resistive dissipation, respectively. The generation due to fluid compression, · , is very weak, and we ignore it throughout. Figure 17 shows profiles of the terms in Equation (15) averaged over a short time interval during fourfold-wreath cycling. The time interval chosen is not unique, and the profiles in Figure 17 are representative of the generation terms at any time the fourfold-wreath cycle is dominant. The governing term is clearly the mean shear, or Ω effect ( f MS ). The fluctuating shear ( f FS ) is the weakest term overall and has no obvious pattern in its spatial distribution. The total advection ( f TA ) is weak and also has no cellular structure (compare to Figure 1(d)), suggesting that transport of B f by the meridional circulation is negligible. The resistive dissipation ( f RD ) plays some role (especially near the inner boundary) and has a sign structure opposite to that of á ñ f B , acting to destroy whatever field is already present.
The mean-shear production of á ñ f B has a similar pattern to á ñ f B itself, but with the four wreaths displaced closer to the equator, which tends to drive the á ñ f B wreath structure equatorward. Furthermore, like the fourfold wreaths themselves, the wreaths of f MS reappear at midlatitudes after they disappear near the equator. In Figure 17(a), for example, the negativef MS wreath straddling the (- 15 ) latitude in the southern hemisphere is associated with a weak negativef MS wreath forming at ∼(- 40 ) latitude. Positive B f thus forms at high latitudes where originally there was negative B f , creating a polarity reversal.
It is unclear exactly what stops the equatorward migrationi.e., why the wreaths have significantly reduced amplitude in the   10 latitude band straddling the equator. In the Sun, there is presumably only one wreath per hemisphere, each with opposite polarity (as inferred from the solar butterfly diagram). When the two wreaths meet at the equator, they annihilate, as in magnetic reconnection. The same process cannot occur in the fourfold-wreath cycle because the wreaths are often equatorially symmetric.
We now investigate the origins of the displacement of the shear-production pattern relative to the azimuthal magnetic  To investigate the source of the meridional field B m , we turn to the radial component of the longitude-averaged Equation (5), which we write as Here, "MT" and "FT" refer to the mean transport (magnetic field generation from the mean electromotive force-or emf -á ñ´á ñ v B ) and the fluctuating transport (or generation from the fluctuating emf á ¢´¢ñ v B ), respectively. Figure 19 shows the dominant generation terms for á ñ B r alongside á ñ B r itself. We only show the fluctuating transport FT r , which is an order of magnitude larger than both MT r and RD r . Similar to the generation of á ñ f B shown in Figure 6, the The fluctuating transport is similar in spirit to the α effect from mean-field theory, which represents the emf caused by the helical action of convection on the magnetic field. In this respect, the fourfold-wreath cycle is characteristic of an α-Ω dynamo. However, our simulations are distinct from most mean-field dynamo models in that the generation of both meridional field and azimuthal field occurs in the same location in the CZ, with the meridional circulation playing essentially no role. This stands in contrast to the F-T paradigm, in which the α effect creates meridional field at high latitudes that is then pumped deep into the CZ by the meridional circulation to seed the Ω effect in the following cycle.

Dynamics of the Partial-wreath Cycle
The physical mechanisms responsible for the system transitioning between the fourfold-wreath and partial-wreath states remain unclear. Nevertheless, we can understand the maintenance of the dynamo during the partial-wreath state in terms of the dominant induction terms. Furthermore, even though the partial-wreath structure is inherently nonaxisymmetric, it is sufficient to perform a mean-field analysis-i.e., identify the dominant longitude-averaged induction terms of Equations (15) and (19). Because the longitude-averaged field á ñ f B has the same polarity as the dominant partial wreath, there is a one-to-one correspondence between cycles in á ñ f B and partial-wreath cycles.
The mean-field generation terms during partial-wreath cycling vary rapidly in space and time. Therefore, significant spatial and temporal averaging must be performed to extract a meaningful signal. Figure 20 shows the spatially averaged contributions to the production of á ñ f B in case D3-1 from the mean shear (which is clearly dominant) and the other terms, as well as á ñ f B itself, during the interval of partial-wreath cycling depicted in Figure 10. All three profiles are roughly sinusoidal. The profile of mean-shear production slightly leads the profile of á ñ f B in time, indicating that the rise and fall of á ñ f B is directly caused by the mean-shear production. The contribution from the other induction terms (which is dominated by the resistive dissipation) is ∼180°out of phase with á ñ f B , meaning the magnetic diffusion attempts to destroy whatever field is already there.
The production of á ñ B m (which in turn affects the production of á ñ f B through the mean shear) is highly chaotic in time during the partial-wreath cycle, and so instead of plotting the production terms in Equation (19)  at mid-depth, where the mean is taken over the time interval (3000, 3300) P rot and the latitude band between 0°and 20°s outh. Similar relationships hold for the production of á ñ q B , which is also dominated by the fluctuating emf.

Bistability in the Context of Solar Observations
Continuous, full-disk observations of the photospheric magnetic field, produced by MDI/SOHO and HMI/SDO, have greatly constrained certain features of the solar dynamo and yielded a number of surprises. Stenflo & Kosovichev (2012) analyzed a large set of MDI magnetograms from 1995 to 2011 and found that Joy's law (which states that emerging sunspot pairs are tilted on average such that the leading spot is closer to the equator, and that the tilt is proportional to heliographic latitude) holds for active regions both large and small. Furthermore, there was no observed dependence of tilt angle on the total active-region flux (or region area), implying that the tilts are inherent to buoyant flux ropes in the solar interior. This stands in contrast to the argument that the tilts are established by the Coriolis force as the flux ropes rise through the CZ (D'Silva & Choudhuri 1993;Fan et al. 1994).
The view that the interior flux ropes are associated with a preferred sense of tilt is interesting in light of the fourfoldwreath cycle. Each of the fourfold wreaths possesses a particular sense of twist as in Figure 18(a). If loops were to spontaneously break off from the wreaths (this does not occur in our simulations because the magnetic diffusivity is too high), they would be endowed with a tilt from the twisting magnetic field lines in the wreath. Nelson et al. (2014) achieved the separation of magnetic loops from interior wreaths and found that the sample of loops reaching the surface made contact with Joy's law. This supports the idea that a fourfold-wreath-like structure in the Sun could establish systematic tilts associated with flux ropes in the deep interior. Stenflo & Kosovichev (2012) further reported that a statistically significant fraction of active regions (5%-25%, depending on region size) violate Hale's polarity law, strongly suggesting that there must be multiple reservoirs of interior magnetism from which buoyant flux ropes originate. Additionally, the "anti-Hale" active regions could appear in the same latitude band as "Hale" active regions, suggesting close Note. "KE" and "ME" refer to the kinetic and magnetic energy densities, respectively. Units are cgs, with a common exponent divided out. The subscripts f, m, and c refer to the energy contributions from the mean azimuthal field components (á ñ f v and á ñ f B ), the mean meridional field components (á ñ v m and á ñ B m ), and the convective components ( ¢ v and ¢ B ), respectively. The fraction of the total energy (no subscripts; magnetic and kinetic energies considered separately) from each contribution is shown as the percentage in parentheses. proximity in latitude for the reservoirs. This result agrees with the geometry of both the fourfold-wreath structure (whose wreaths come in opposite-polarity pairs with wreaths in a pair separated in latitude by ∼15°) and the partial-wreath structure (with opposite-polarity partial wreaths at approximately the same latitude).
Li (2018) performed another extensive analysis of sunspot groups from 1996 to 2018 using both the MDI and HMI magnetograms, also reporting that a sizable fraction (∼8%) of sunspot pairs violate Hale's polarity law and confirming the result obtained by Stenflo & Kosovichev (2012). In addition, Li (2018) found that the anti-Hale sunspots had different properties from Hale sunspots, namely smaller bipole separation and weaker total magnetic flux. These results further support the idea that the Hale and anti-Hale sunspots may come from different interior reservoirs.
It is also possible that the anti-Hale active regions form an essential part of the solar dynamo. Hathaway & Upton (2016) find that modeling the flux transport of active regions observed by SDO during the current solar cycle can be used to predict the amplitude and hemispheric asymmetry of the following cycle. Central to such predictions is the initial distribution of activeregion tilts, which would be significantly modified by the presence of anti-Hale active regions. Although it is unlikely that the real interior solar magnetic field matches the fourfoldwreath and partial-wreath structures in detail, such configurations offer two ways in which violations of Hale's polarity law could occur.
We now discuss the asymmetric features of solar activity, many of which are captured by the partial-wreath cycle. For well over a century, sunspots have been observed to emerge preferentially at particular solar longitudes (e.g., Svalgaard & Wilcox 1975;Bogart 1982;Henney & Harvey 2002;Maunder 2016). This phenomenon, referred to as active longitudes, remains poorly understood.
It is common for active longitudes to come in pairs, separated by about  180 in longitude (e.g., Bai 2003). The pairs are also usually found in a single hemisphere and are part of the global hemispheric asymmetry in magnetic activity. Finally, the pairs are associated with opposite polarities of the azimuthal magnetic field on opposite sides of the Sun, as revealed by synoptic maps from the Wilcox Solar Observatory (Mordinov & Kitchatinov 2004). Thus, it is possible that the pairs are associated with nests of Hale and anti-Hale active regions on opposite sides of the Sun, although a systematic study of the MDI/HMI magnetograms would be needed to confirm this result.   Active longitudes and their associated hemispheric asymmetry bear a striking resemblance to the partial-wreath structure in our simulated dynamos. The partial wreaths preferentially occupy one hemisphere at a time. They come in opposite-polarity pairs with the strongest magnetism on opposite sides of the sphere (see the Mollweide projection in Figure 2(b) and the 3D volume renderings in Figure 11). The time-longitude diagrams for the partial wreaths in case D3-1 ( Figure 12) show that each wreath is roughly stationary once an appropriate rotating frame is chosen, similar to the quasi-rigid structure formed by active-longitude pairs in the Carrington reference frame (Berdyugina & Usoskin 2003). Figures 11 and 12 would also indicate cyclic behavior in which the strongest active longitude flips by  180 over the partial-wreath cycle period. Berdyugina & Usoskin (2003) reported such behavior in observations of solar sunspot number, wherein the active longitudes appeared to flip between opposite sides of the Sun every ∼3.7 yr, or every ∼53 solar Carrington rotations. Our simulation suite displays similar properties in the relative energy content of the partial wreaths ( Figure 13), showing that the dominant partial wreath switches with a quasi-regular period of 81.5 P rot .
The hemispheric asymmetry in solar activity can be seen in a wide variety of activity indicators, such as number of flares, sunspot number, sunspot area, and photospheric magnetic flux (Schüssler & Cameron 2018). Many studies have been performed on the periodic behavior of the asymmetry, with most agreeing on a long-term period of ∼50 yr and a short-term period of ∼10 yr (see Hathaway 2015 and references therein). Our simulations each possess a significant hemispheric asymmetry (Figure 14(a)), although there is no obvious associated cycle (Figure 14(b)). Nevertheless, our simulations make contact with observations in that the hemispheric asymmetry occurs as a long-term modulation of another cycle, which in our simulations is that of the partial wreaths.
In summary, the dynamo cases presented here show similarities to several features of solar activity, namely the sunspot cycle, opposite-polarity magnetic regions in close proximity, active longitudes, and hemispheric asymmetry. We do not suggest that all these observed features of the Sun's activity are a result of fourfold-wreath and partial-wreath structures existing in the solar CZ. Rather, our simulations provide insights into to how a solar-like CZ can support dynamos exhibiting a wide range of cyclic behavior.

Conclusions
We have explored three dynamo simulations of a 3D solarlike CZ at varying magnetic Prandtl numbers. These dynamo cases display clear signs of bistability. In particular, there is a fourfold-wreath structure that exhibits equatorward propagation and regular polarity reversals. There is also a uni-hemispherical partial-wreath pair that is roughly stationary and whose dominant polarity changes cyclically, but with an irregular period. The fourfold-wreath cycle captures aspects of the 22 yr sunspot cycle (namely, the solar butterfly diagram), whereas the partial-wreath cycle is reminiscent of several of the observed features of active longitudes and hemispheric asymmetry.
We have examined the dynamics of the fourfold-wreath cycle in detail and find they are approximately described by an α-Ω dynamo. The mean shear, or Ω effect, is responsible for generating the wreaths. The regular polarity reversals and equatorward propagation are caused by the inherent geometry of the wreaths (namely the meridional field winding counterclockwise for positive B f and clockwise for negative B f ) coupled to a solar-like differential rotation. Generation of the meridional field occurs through the fluctuating emf, which is similar to an α effect, but fully nonlinear. Both effects occur in the same localized region of the CZ (with the meridional circulation playing basically no role), unlike in the F-T dynamo framework. As a whole, the production terms for the fourfoldwreath cycle are similar to those of model K3S in Augustson et al. (2015), although in that work, significant modulation of the differential rotation played a large role, whereas in our dynamo cases, the differential rotation stays mostly constant with time.
The partial-wreath cycle has longitude-averaged dynamics similar to an α-Ω dynamo as well, but it is not fully clear how the system transitions between the regularly cycling fourfold wreaths and the uni-hemispherical partial-wreath pair. It is possible that the transition is governed by the fourfold-wreath structure changing between equatorially symmetric and antisymmetric states, but more work needs to be done to assess in detail how such a mechanism operates. Equatorial parity also played a key role in model K3S of Augustson et al. (2015), wherein atypical dominance of the radial magnetic field by the even-ℓmodes acted as a precursor to the cycling dynamo entering a "grand minimum." Similarly, in the global convective dynamo simulation of Raynaud & Tobias (2016), the interaction of hydromagnetic modes with different symmetries led to transitions between distinct cycling Figure 20. The dominant generation terms appearing in Equation (19), as well as á ñ f B , when the partial-wreath structure is cycling in the south, during the same interval shown in Figure 10). Quantities are averaged over radius and between 0°and 20°south in latitude. The generation terms are broken up into mean shear and "other" (i.e., ¶á ñ ¶f f B t MS ) and multiplied by 5 P rot to put them on the same scale as á ñ f B . The contribution to "other" is dominated by resistive dissipation. dynamo states. The magnetic Prandtl number may also be significant in determining whether the system is unstable to the partial-wreath state; Nelson et al. (2013) found magnetic structures resembling uni-hemispherical partial wreaths in their high-Pr m simulations but not low Pr m . Our work suggests that spherical-shell convection can support complex dynamos capable of exhibiting vastly different global structures for the magnetism simultaneously. In particular, the bistable properties of our dynamo cases are consistent with a solar dynamo containing two cycles at once: one associated with the observed butterfly diagram and one associated with active longitudes and hemispheric asymmetry. We recognize that our simulations here are just starting to explore the complexities that might be at work in the highly nonlinear and turbulent dynamics of the Sun. The range of scales that we can capture and the diffusivities we use are still far removed from the solar regime. However, we have found that some reasonable large-scale behavior emerges from our dynamo models, and these may broaden our insights into the physical processes occurring deep in the interior of the Sun. for cases D3-1, D3-2, and D3-4, respectively.