© The Authors 2025

1. Introduction

Widespread solar energetic particle (SEP) events, spanning more than 180° around the Sun, remain poorly understood in heliophysics. They are often associated with large coronal mass ejections (CMEs) and the shocks they generate (Rouillard et al. 2012; Lario et al. 2014), but the exact mechanisms that enable particles to spread so widely are still debated (e.g. Kollhoff et al. 2021; Rodríguez-García et al. 2021; Kouloumvakos et al. 2022; Dresing et al. 2023). One possibility involves wide, possibly circumsolar, shock waves (e.g. Hou et al. 2022), which can accelerate particles throughout extensive regions of the heliosphere (Gómez-Herrero et al. 2015; Lario et al. 2016). Observations of large shock extents (e.g. Nitta et al. 2013; Kwon & Vourlidas 2017) support this view. Alternatively, SEPs may originate from localized regions and spread via cross-field transport in the corona and interplanetary (IP) medium (e.g. Dresing et al. 2012; Strauss et al. 2017).

Here, we focus on the role of wide shocks in driving widespread SEP events. A key question is how far and for how long shocks – which are often significantly wider than their parent CMEs – can persist. Although the solar wind and corona are complex, shocks typically evolve as they propagate, with some regions remaining CME-driven while others expand beyond direct CME influence (e.g. Corona-Romero & Gonzalez-Esparza 2011; Warmuth 2015). This transition begins at the shock flanks, which decouple first, while the central region remains CME-driven. Over time, as the CME slows and its ability to drive the shock diminishes, even the central part of the shock can detach, leaving behind a fully decoupled wavefront (e.g. Pinter & Dryer 1990).

Observations by Kwon & Vourlidas (2018) highlight that shock strength can remain substantial well beyond the region directly driven by the CME apex, indicating that wide CME-driven shocks can accelerate energetic particles over larger spatial and temporal scales than the CME width alone would suggest. Shock modelling and reconstructions generally support these observations (e.g. Kouloumvakos et al. 2019; Jebaraj et al. 2021), though there are exceptions (e.g. Lario et al. 2017; Kollhoff et al. 2021), illustrating the complexity of shock evolution and SEP transport.

A related issue is the fate of shocks as they propagate into IP space. Once partially or fully decoupled from their drivers, shocks rely on the energy and momentum they gained during their driving phase (e.g. Corona-Romero & Gonzalez-Esparza 2012). Their subsequent evolution is reminiscent of a blast wave (e.g. Smart & Shea 1985; Pinter & Dryer 1990; Corona-Romero & Gonzalez-Esparza 2011; Li et al. 2024) and depends strongly on the upstream conditions, such as the solar wind density, magnetic field structure, and speed. In some cases, shocks can maintain sufficient strength to accelerate particles over large distances and for a long time, while others weaken relatively quickly (Liu et al. 2017).

A key factor in this evolution is the decrease in the local fast magnetosonic speed with heliocentric distance. For instance, Echer (2019) show that IP shocks near Jupiter’s orbit can exhibit higher Mach numbers on average than those at 1 au, not because of additional energy from the CME, but because the magnetosonic speed at 5 au is only about 60% of its value at 1 au. Consequently, shocks that appear modest closer to the Sun can seem stronger farther out simply due to the changing local solar wind environment. Understanding these factors is essential for accurately modelling how long shocks can persist and the conditions under which they can continue to accelerate SEPs throughout the heliosphere.

In this work we explored the possibility that the 13 March 2023 event involved two nearly simultaneous CMEs whose combined eruption drove a wide IP shock. We examine whether the flanks of this shock could have decoupled from the ejecta and evolved into freely propagating waves as the structure expanded through the heliosphere. To investigate this scenario, we used the European Heliospheric Forecasting Information Asset (EUHFORIA; Pomoell & Poedts 2018), which solves 3D, time-dependent ideal magnetohydrodynamic (MHD) equations from 0.1 au onwards. As part of this study, we upgraded EUHFORIA to enable the direct injection of shocks, rather than relying solely on ejecta-driven shock waves (see Appendix A). In our simulations, a magnetic ejecta continues to drive the shock at its apex, while the shock flanks can propagate freely. Under these modelling assumptions, the shock can maintain its strength out to large radial distances, partly due to the decreasing fast magnetosonic speed farther from the Sun. Our results show promising agreement with in situ observations at various longitudes and radial distances, suggesting that such wide shocks can persist well beyond the CME core.

2. The 13 March 2023 event: Observations

On 13 March 2023 at approximately 03:15 UT, significant solar activity on the Sun’s far side produced a fast and wide shock wave. Parker Solar Probe (Fox et al. 2016), positioned closest to the eruption site, observed an extraordinary shock wave followed by a magnetic cloud at 0.24 au. Jebaraj et al. (2024) reported that this shock was among the fastest in situ shocks ever recorded (V_sh ∼ 3000 km s⁻¹).

The eruption’s source region remains uncertain due to limited far-side imaging. However, Dresing et al. (2025) identified two candidate active regions (ARs) by analysing full-disk observations before and after the event. One AR (AR1) was located in the northern hemisphere, initially designated AR13229 during Carrington rotation 2267 and reappearing as AR13258 in the subsequent rotation. The other AR (AR2), in the southern hemisphere, was originally AR13236 and later evolved into a group of ARs (AR13256, AR13257, and AR13259) during Carrington rotation 2268. Stonyhurst coordinates place these ARs at approximately −109° (AR1) and −169° (AR2) longitude at the time of eruption. The estimated positions of these ARs are marked in Fig. 1, along with the locations of spacecraft at the time of the event. For a detailed analysis of AR identification and evolution, see Dresing et al. (2025).

Fig. 1. Spacecraft constellation in the solar equatorial plane, showing nominal Parker spiral field lines connecting each spacecraft to the Sun. The field lines were derived using measured solar wind speeds for STEREO-A, Earth, and Parker Solar Probe, while for Solar Orbiter and BepiColombo, where no solar wind plasma data were available at the time of the event, solar wind speeds from the EUHFORIA simulation are used. The two arrows indicate the estimated locations of the ARs from which the CMEs erupted. The dashed pink line represents the combined extent of the two shock waves injected in the EUHFORIA simulation. This figure was made using the Solar-Mach tool (Gieseler et al. 2023).

Figure 2 presents coronagraph observations of the CME and the CME-driven shock wave in white light, taken by the Large Angle and Spectrometric Coronagraph Experiment (LASCO; Brueckner et al. 1995) C2 camera on board the SOlar and Heliospheric Observatory (SOHO; Domingo et al. 1995) and the COR2 coronagraph (Howard et al. 2008) on board the Solar TErrestrial RElations Observatory Ahead (STEREO-A; Kaiser et al. 2008) spacecraft. The event likely involved two simultaneous CMEs on the Sun’s far side – one moving north and the other south (see panels a and b) – originating from AR1 and AR2, respectively. Both CMEs were enveloped by a bright white-light front identified as the outer shock wave. This shock is clearly visible in the coronagraph images, and the deflection of several streamers on the west limb (see S1D and S2D in panels d–f) confirms a wide disturbance propagating in the corona.

Fig. 2. Selected snapshots of the CME evolution on 13 March 2023, as seen in coronagraphic observations from LASCO and STEREO-A. Each panel shows an overlay of a plain coronagraphic image and the corresponding running difference image at the same time from LASCO-C2 (panels a and b) and COR2-A (panels c to f). At selected frames we note the CMEs and outline the location of the CME-driven shock wave. The interaction of the shock with streamers (S1 and S2) are also noted in panels d to f (streamers deflections S1D and S2D). The various features are discussed in detail in Sect. 2.

The shock was also linked to an intense, circumsolar SEP event (see Dresing et al. 2025), with particle intensity increases recorded at Parker Solar Probe, the first Sun–Earth Lagrange point (L1), Solar Orbiter (Müller et al. 2020), STEREO-A, BepiColombo (Benkhoff et al. 2021), and possibly even Mars. Notably, Earth’s nominal Parker spiral connection lies about 120° in longitude from AR2 (Fig. 1), yet it still experienced a rapid SEP onset. According to Dresing et al. (2025), the observed long-lasting anisotropy at L1 suggests either very efficient particle transport across large longitudinal separations close to the Sun or a direct magnetic connection to a wide shock wave, likely the eastern flank of the CME erupting from AR2. Furthermore, each spacecraft eventually recorded an energetic storm particle (ESP) event when a shock arrived, despite the broad range of heliolongitudes involved. As shown in the sections that follow, a possible explanation is that the western flank of the shock wave originating from AR1 propagated to the far-side spacecraft, suggesting that the IP shock associated with the 13 March 2023 event was not only exceptionally fast but also strikingly wide.

3. EUHFORIA setup for the 13 March 2023 event

To investigate the potential role of a wide shock wave in the 13 March 2023 eruption, we performed a simulation using the EUHFORIA model (Pomoell & Poedts 2018; Poedts et al. 2020). EUHFORIA consists of two primary components: a coronal model that derives solar wind and magnetic field conditions at 0.1 au from photospheric magnetograms, and a heliospheric model that propagates the plasma and magnetic structures to larger radial distances by solving the ideal MHD equations. In this work, we employed the Air Force Data Assimilative Photospheric Flux Transport (ADAPT; Arge et al. 2010; Hickmann et al. 2015) model in conjunction with Global Oscillation Network Group (GONG; Harvey et al. 1996) magnetograms to specify the coronal conditions. Specifically, we used the first GONG-ADAPT realization produced at 00:00 UT on 13 March 2023, and then set the MHD boundary conditions at 0.1 au following Pomoell & Poedts (2018), with two modifications: (i) Increasing the reference magnetic field at 0.1 au from 300 nT to B_fsw = 500 nT, to achieve better agreement with in situ spacecraft data. (ii) Raising the minimum solar wind speed from 275 to 350 km s⁻¹, ensuring that the solar wind remains superfast at the inner boundary (0.1 au), as required by the MHD module of EUHFORIA (e.g. Samara et al. 2024). Our computational domain extends to 3 au, with a radial grid spacing of Δr = 0.6 R_s and an angular resolution of 1° in both latitude and longitude.

To better capture the heliospheric conditions prior to the eruption of the primary CMEs on 13 March 2023, our EUHFORIA simulation includes five pre-event CMEs and one post-event CME. These CMEs were selected based on their potential impact on the heliospheric environment surrounding the main event. Specifically, we included CMEs that erupted within three days before the main event and met at least one of the following criteria: (i) Earth-directed CMEs with projected speeds exceeding 350 km s⁻¹, or (ii) CMEs with angular widths greater than 90° and speeds above 500 km s⁻¹, as listed in the CDAW SOHO LASCO CME catalogue¹ (Yashiro et al. 2004). The kinematic parameters of these CMEs are detailed in Table 1. To better constrain these parameters, we applied the graduated cylindrical shell model (Thernisien et al. 2006; Thernisien 2011) to multi-point coronagraph observations from STEREO-A and SOHO, reconstructing the CME properties from approximately 3.5 to 20 R_⊙.

To model the 13 March 2023 eruption, we injected two spheromak CMEs at 0.1 au at the inferred locations of AR1 and AR2 (see Fig. 1). Each spheromak was initialized with a radius of 12 R_⊙, corresponding to a ∼30° half-width. Prior to the spheromak injections, we introduced two wide shock waves at the inner boundary by applying Rankine–Hugoniot (RH) jump conditions, as described in Appendix A. These shocks were assigned a half-width of 110° and a uniform speed of 3000 km s⁻¹, consistent with estimates by Jebaraj et al. (2024). The central region of the wave fronts remains CME-driven, while the flanks propagate freely. The shock centres, like in the spheromaks, coincide with the estimated AR locations and are placed 60° apart, resulting in significant overlap between their wave fronts. In these overlapping regions, the shocks are treated as a single wave, determined by the shock speed and local upstream solar wind conditions. This setup aligns with coronagraph observations, which suggest a shock envelope encompassing both CMEs (Fig. 2). The total azimuthal extent of the injected shock waves is 280°, as shown in Fig. 1.

We adopted a 110° half-width for the shock from AR1, ensuring that the western tip of the shock flank travels in the direction of L1. Applying the same half-width to the shock from AR2 places its eastern flank at 81° in azimuthal angle (Stonyhurst coordinates), which remains approximately 10° away from L1’s nominal magnetic connection to the solar surface (see Dresing et al. 2025, or Fig. 1). Given the observed anisotropic SEP event onset at L1, this suggests that the shock from AR2 may have been even wider than assumed. Moreover, while we impose a constant shock speed of 3000 km s⁻¹ across the entire wave front for simplicity, the actual speed is likely lower at the flanks than at the CME-driven centre. Nevertheless, this setup allows us to investigate the evolution of wide shock flanks in a preconditioned heliosphere, gain further insights into the shock’s large-scale structure, and make comparisons with observations from multiple spacecraft.

4. Simulation results

Figure 3 presents a comparison between the EUHFORIA simulation results and observations from various spacecraft, including Parker Solar Probe, Solar Orbiter, STEREO-A, and spacecraft at L1. This figure displays a range of in situ measurements alongside the corresponding quantities from the EUHFORIA simulation. From top to bottom, each panel shows the solar wind proton speed, the proton number density, the thermal pressure, the magnetic field magnitude, the Geocentric Solar Ecliptic (GSE) latitudinal magnetic field angle, and the GSE azimuthal magnetic field angle. The arrival of the shock wave is clearly identifiable at each spacecraft by a sudden jump in most of these panels and is marked by a dotted line. Comparing these in situ observations with the EUHFORIA simulation results, we find a remarkable agreement. Both the arrival time and the amplitude of the shock wave align well across all spacecraft.

Fig. 3. Comparison of observed (blue) and simulated (orange) solar wind parameters at various spacecraft locations: Parker Solar Probe (top left), Solar Orbiter (top right), L1 (ACE and WIND; bottom left), and STEREO-A (bottom right). For each spacecraft, the panels display: (a) solar wind speed, (b) proton number density, (c) thermal proton pressure, (d) magnetic field magnitude, (e) magnetic latitudinal angle, and (f) magnetic azimuthal angle.

An important observation from the multi-spacecraft analysis of in situ data is that, apart from Parker Solar Probe, no other spacecraft detected clear signatures of a flux rope downstream of the shock. This not only aligns well with our simulation results but also supports the hypothesis of a shock flank encounter without any CME signatures to follow. A typical characteristic of a flux rope is a smooth and coherent rotation in the magnetic field angles, yet no such structure is evident in the in situ magnetic field data. Solar Orbiter is the only spacecraft to exhibit some rotation in the magnetic field at the end of the enhanced field region; however, this coincides with an increase in plasma density, which is not a typical flux-rope signature. Furthermore, the enhanced magnetic field strength observed downstream of the shock at multiple spacecraft, along with fluctuations in the magnetic field direction, is expected in a turbulent shock sheath and does not necessarily indicate the presence of a flux rope.

This is further supported by the EUHFORIA simulation, which, for example, shows variations in the magnetic field angles at STEREO-A, despite no magnetic cloud passing the spacecraft in the simulation. Moreover, the fact that the simulated magnetic field agrees relatively well with observations – both in magnitude and duration – at spacecraft encountering the modelled shock flank further suggests that no flux rope crossed these spacecraft. This strengthens the interpretation that the observed magnetic field variations are a consequence of the turbulent shock sheath rather than the passage of an underlying magnetic ejecta structure. Alternatively, these variations could result from glancing encounters with the less-organized flanks of the ejecta, where the magnetic field structure is more complex and lacks the coherent rotation expected for a magnetic cloud (e.g. Rodríguez-García et al. 2022).

Figure 4 shows a snapshot from the EUHFORIA simulation on 14 March at 05:57 UT, illustrating the solar wind speed (top row) and the divergence of solar wind velocity, ∇ ⋅ V (bottom row). The left column displays the solar equatorial plane, while the right column shows a meridional slice passing through Earth. In the solar equatorial plane, the shock wave is distinctly visible as a pronounced jump in speed, encircling the plot’s centre. Two spheromaks representing the main event appear as extended regions of enhanced velocity, concentrated primarily in the third quadrant. In this study, we do not focus on these spheromak CMEs but instead analyse the freely propagating wave, also visible in the other three quadrants of the EUHFORIA simulation.

Fig. 4. Snapshot of the EUHFORIA simulation on 14 March at 05:57 UT. Left panels: Solar equatorial plane. Right panels: Meridional slices containing Earth. Top row: Solar wind speed. Bottom row: Divergence of the solar wind velocity ∇ ⋅ Vsw. Dashed lines indicate projections of magnetic field lines connecting to various spacecraft. A video version of this figure is available online.

The bottom row of Fig. 4 displays ∇ ⋅ V and highlights the main-event shock wave. Here, ∇ ⋅ V < 0 represents converging flows, typical of compression regions or shock fronts, while ∇ ⋅ V > 0 indicates diverging flows, characteristic of the expanding solar wind. Downstream of the main-event shock, the plasma is characterized by ∇ ⋅ V ≫ 0, indicating a strongly expanding flow.

At the time of the snapshot shown in Fig. 4, the main-event shock had already overtaken two earlier CMEs, launched on 12 and 13 March, respectively (see Table 1). Initially, the shock spanned 280° in longitude, from 1° to −281° ( = 81°). However, due to the lateral expansion of its flanks, the disturbance has since become quasi-circumsolar in the equatorial plane, with both the eastern and western flanks expanding by over ∼40° and thus merging around 40° longitude. Additionally, our simulation likely underestimates the lateral expansion of the shock, as it is injected with a purely radial velocity at 0.1 au. As a result, a narrower initial shock half-width, such as 70° instead of 110° from AR1, should, in principle, still have allowed the shock to reach L1.

In addition to the main-event shock, several other shock waves are visible, associated with preceding CMEs included in the simulation (see Table 1). For example, the meridional slice in Fig. 4 shows a pre-event CME whose northern flank is expected to pass Earth around 06:00 UT on 14 March (see also Fig. 3). In the solar equatorial plane, the three earliest pre-event CMEs have merged into a single, broad compression front. This type of merging is reminiscent of what likely occurred during the main event itself, when two fast CMEs erupted almost simultaneously, producing the exceptionally wide shock wave analysed in this work.

Finally, the post-event CME, propagating through the shocked solar wind behind the main-event shock, is also visible in Fig. 4. As shown more clearly in the video version of the figure, this post-event shock initially catches up with the main-event CME, but ultimately dissipates before ever overtaking the main-event shock.

Focusing on the circumsolar wave, three key characteristics stand out from this 3D analysis: (i) the freely propagating flanks travel significantly slower than the central, spheromak-driven portion of the shock, (ii) the wave front exhibits considerable non-uniformity and deformation along its extent, and (iii) it is predominantly characterized by ∇ ⋅ V ≪ 0, except at the very tips of the flanks, suggesting that it likely persists as a shock wave rather than transitioning into a simple compression wave.

To confirm that the disturbance is a shock wave, we employ the shock wave tracer described in Wijsen et al. (2022) and calculate its properties such as the fast-magnetosonic Mach number, M_fms. The top panel of Fig. 5 shows the shock surface on 14 March at 05:58 UT, colour coded according to M_fms. As indicated, M_fms > 1 throughout, verifying that the wave is indeed a shock. A 360° video of the shock is also available online. As expected, the video reveals that the shock is strongest – i.e. characterized by the highest M_fms – where the spheromaks are driving the wave. However, even at the flanks, where the wave is not driven by the spheromaks, the wave remains a shock until the end of the simulation, at which point the shock has travelled beyond 2 au.

Fig. 5. Main-event shock surface, along with the magnetic field lines connecting Earth (green), STEREO-A (red), Solar Orbiter (blue), BepiColombo (orange), and Parker Solar Probe (purple). In the top panel, the shock surface is colour-coded by the fast magnetosonic Mach number (Mfms), while in the bottom panel, it is colour-coded by the shock obliquity angle (θBn). The associated online video provides a 360° view around the shock surface.

The pronounced deformation of the shock results from the non-uniform upstream solar wind, a phenomenon previously discussed for CME-driven shocks in, for example, Wijsen et al. (2023). This deformation gives rise to a shock with highly variable geometry, as demonstrated in the bottom panel of Fig. 5. This panel shows the shock obliquity θ_Bn, defined as the angle between the shock normal and the upstream magnetic field. Shock obliquity is believed to play a crucial role in the acceleration of energetic particles. For example, diffusive shock acceleration is thought to be most efficient for quasi-parallel shocks, where the obliquity angle is near 0° (e.g. Malkov & Drury 2001, and references therein).

To further illustrate the evolution of the freely propagating shock wave, Fig. 6 shows its progression along two radial directions in the solar equatorial plane. The first direction (blue symbols), at a longitude of 0° degrees, aligns with the x-axis in Fig. 4. The second direction (green symbols), at 135° longitude, corresponds to the radial line in the second (upper left) quadrant of the solar equatorial plane in Fig. 4. The 135° direction was selected because no pre-event CME propagates along this path, meaning that the shock interacts only with the unperturbed solar wind.

Fig. 6. Evolution of the main-event shock flanks along two radial directions in the solar equatorial plane. The blue symbols represent the line at 0°, aligned with the positive x-axis in Fig. 4, while the green symbols represent a line rotated 135° anti-clockwise from the positive x-axis. The panels present: (a) hours elapsed since the blast wave injection, including power-law fits and observed spacecraft shock crossing times; (b) shock speed; (c) upstream solar wind speed; (d) upstream solar wind density, scaled with r2; (e) upstream fast-magnetosonic speed; and (f) fast-magnetosonic Mach number (Mfms).

Figure 6a presents the time elapsed since the initial blast wave injection on 13 March at 05:00 UT as a function of radial distance from the Sun, with observed shock crossing times for various spacecraft also indicated. These observations align well with the simulation, reinforcing the accuracy of the model. The panel also includes power-law fits of the form t ∝ r^α. Near the model’s inner boundary (0.1–0.4 au), we find power-law exponents of r^{2.358 ± 0.19} and r^{2.29 ± 0.17} for the 0° and 135° trajectories, respectively. Beyond 0.4 au, the expansion slows, with exponents of r^{1.46 ± 0.01} and r^{1.57 ± 0.01} for the same respective trajectories. These values are consistent with the analytical behaviour of freely propagating pressure waves. In an idealized scenario where energy is deposited instantaneously (e.g. an explosion), the shock initially decays rapidly before transitioning into a regular pressure wave. This problem, well-studied in uniform hydrodynamic media, is described by the Sedov–von Neumann–Taylor (SVT) solution (Sedov 1946; Taylor 1950b,a; von Neumann 1963), which predicts t ∝ r^3/2 for a medium where density decreases as 1/r² (see also Cavaliere & Messina 1976 and our Appendix B). Similar dependences for freely propagating shock stages in the solar wind have been identified in previous studies. Smart & Shea (1985) examined the evolution of IP shocks and found that once a CME-driven shock ceases to be actively driven by the ejecta, its propagation follows a power-law decay akin to the SVT blast wave solution. Pinter & Dryer (1990) further analysed decoupled IP shocks and demonstrated that their radial evolution can be described by self-similar solutions, with a transition from a piston-driven phase to a freely propagating state. However, the SVT solution assumes a uniform, expanding upstream plasma with negligible pressure compared to the downstream, so an exact match is not expected. Notably, near the inner boundary, particularly along the 0° trajectory, the slope is steeper due to interaction with a pre-event CME.

Panels (c) and (d) of Fig. 6 show the variability in the upstream solar wind speed and density, respectively. Along the 135° line, a significant density enhancement occurs around 0.75 au, corresponding to stream interaction regions encountered by the shock wave as it propagates outwards. Along the 0° radial line, two prominent density peaks are evident near 0.2 au and 1.25 au, associated with two pre-event CMEs overtaken by the freely propagating shock.

The evolution of the fast magnetosonic speed, V_fms, and Mach number, M_fms, is shown in panels (e) and (f), respectively. Notably, V_fms decreases almost monotonically with increasing heliocentric distance, a trend that plays a crucial role in sustaining the shock wave as it propagates outwards. Since solar wind plasma that was initially shocked at smaller radial distances encounters progressively lower magnetosonic speeds, it can continue to reinforce the shock structure at larger distances without requiring additional energy input. As a result, shocks that appear relatively weak near the Sun can seem stronger farther out purely due to the evolving solar wind environment (see also Echer 2019).

This effect is particularly pronounced along the 135° trajectory, where M_fms increases significantly around 0.8 au and 1.6 au. The first increase is driven by the sharp decrease in V_fms as the shock encounters a stream interaction region. The second increase is resulting from the shock propagating into a low-density rarefaction region trailing a fast wind stream. Finally, as shown in panel (f), M_fms consistently exceeds unity along both radial lines, confirming that the disturbance remains a shock wave throughout its propagation from the upper corona to 2 au.

5. Conclusions

In this study we investigated how the flanks of a wide CME-driven shock can persist in the heliosphere as freely propagating waves, using the EUHFORIA MHD solar wind model. We modified EUHFORIA to allow the direct injection of shock waves with specified speeds at the inner boundary by ensuring that their parameters meet the RH jump conditions for the given shock speed. Our key findings are summarized as follows:

1. Sustained, freely propagating shock flanks: EUHFORIA simulations show that the flanks of wide CME-driven shocks can persist beyond 2 au, even in the absence of a continuous CME driving at those flanks. A key factor in this persistence is the decreasing fast magnetosonic speed with heliocentric distance, allowing initially shocked solar wind to sustain the shock farther out without additional energy input.

2. Application to the 13 March 2023 event: We modelled two nearly simultaneous CMEs that drive a single, wide shock. In this configuration, the central shock region near the apex remains CME-driven, whereas the flanks propagate freely. The resulting global shock front, initially injected with a 280° angular width, extends into a quasi-circumsolar wave. The modelled shock arrival times and amplitudes of key plasma parameters (e.g. speed and density) closely match in situ observations from multiple spacecraft spanning various radial distances and longitudes. The 13 March 2023 event represents an extreme case, with shock speeds reaching up to 3000 km s⁻¹ at Parker Solar Probe, leading Jebaraj et al. (2024) to term it an ‘almost astrophysical’ shock. This exceptionally high speed may explain why the shock persisted at such large heliocentric distances. For weaker shocks, the lack of a sustained CME driving at the flanks may lead to a faster dissipation.

3. Implications for widespread SEP events: The persistence of shock flanks may contribute to the spatial extent and duration of widespread SEP events. For extreme cases, like the event discussed in this work, the shock likely needs to propagate well into IP space to help explain ESP events observed across a wide longitudinal range Dresing et al. (2025). Variations in observed particle profiles may thus be influenced not only by transport processes but also by differences in shock geometry and local plasma conditions.

While our simulation results suggest that the 13 March 2023 event can be explained by a single, wide shock, this is only one possible interpretation of the observations. Alternatively, the event could also result from a shock-flank encounter with a preceding CME (see also Dresing et al. 2025). In the presented EUHFORIA simulation, these earlier CMEs either miss the spacecraft or arrive too late (more than 12 hours after the observed shock).

However, heliospheric MHD models such as EUHFORIA and ENLIL (Odstrcil 2003) struggle to accurately reproduce shock-flank encounters (e.g. Bain et al. 2016; Scolini et al. 2019). This is partly because they start their simulations at 0.1 au, underestimating the true shock extent by injecting only the CME ejecta while neglecting the shock wave and sheath that form at lower altitudes. As a result, these models also face difficulties in predicting SEP onset times when an observer connects to the shock flank (see e.g. Wijsen et al. 2022).

By explicitly injecting a shock ahead of the CME ejecta that satisfies the RH jump conditions at the inner boundary, our simulations take a step towards addressing this limitation, providing a more realistic representation of the shock evolution in the heliosphere. In principle, such a shock and sheath should be injected ahead of every sufficiently fast CME. A more comprehensive, though computationally expensive, approach would involve modelling the entire solar eruption from the low corona onwards to fully capture the formation and expansion of the shock (e.g. Manchester et al. 2004; Linan et al. 2025).

Finally, we note that the results obtained in this study rely on an ideal MHD framework, which does not account for certain physical processes that may influence the persistence of freely propagating shock fronts, such as kinetic-scale dissipation. Future work will focus on investigating effects beyond ideal MHD to further improve our understanding of shock evolution and its role in particle acceleration and transport.

Data availability

Movies associated to Figs. 4 & 5 are available at https://www.aanda.org

Acknowledgments

N.W. acknowledges funding from the Research Foundation – Flanders (FWO – Vlaanderen, fellowship no. 1184319N) and from the KU Leuven project 3E241013. Computational resources used for the EUHFORIA simulations presented in this work were provided by the VSC (Flemish Supercomputer Center), funded by FWO – Vlaanderen and the Flemish Government – department EWI. Work in the University of Turku was performed under the umbrella of Finnish Centre of Excellence in Research of Sustainable Space (FORESAIL) funded by the Research Council of Finland (grant No. 352847). I.C.J. and N.D. are grateful for support by the Research Council of Finland (SHOCKSEE, grant No. 346902). We thank the members of the data analysis working group at the Space Research Laboratory of the University of Turku, Finland for useful discussions. We acknowledge funding by the European Union’s Horizon 2020/orizon Europe research and innovation program under grant agreement No. 101004159 (SERPENTINE) and No. 101134999 (SOLER). The paper reflects only the authors’ view and the European Commission is not responsible for any use that may be made of the information it contains. I.C.J. also acknowledges support from the International Space Science Institute (ISSI) in Bern through ISSI International Team project No. 23-575, “Collisionless Shock as a Self-Regulatory System”. A.K. acknowledges financial support from NASA’s NNN06AA01C (80MSFC19F0002 SO-SIS Phase-E and Parker Solar Probe EPI-Lo) contract. E.P. acknowledges support from NASA’s PSP-GI (grant No. 80NSSC22K0349), HGI (grant No. 80NSSC23K0447), LWS (grant No. 80NSSC19K0067), and LWS-SC (grant No. 80NSSC22K0893) programmes, as well as NSF’s PREEVENTS (grant No. ICER-1854790) programme.

References

Appendix A: Modelling freely propagating shock waves

EUHFORIA’s MHD module solves the complete set of ideal MHD equations, including the conservation of mass, momentum, energy, and the induction equation for the magnetic field. The solar wind in the heliosphere is modelled by evolving these equations through a structured mesh. At 0.1 au, the inner boundary, the solar wind velocity is radial and superfast, meaning that the solar wind speed exceeds the local fast magnetosonic speed everywhere. This ensures that all MHD characteristics are directed outwards. Furthermore, CMEs can be introduced at the inner boundary as hydrodynamic plasma blobs — for example, the cone CME or the spheroid model (Scolini & Palmerio 2024) — or as flux-rope-like structures — for example, the spheromak (Verbeke et al. 2019) or the FRi3D (Maharana et al. 2022) descriptions. These CMEs are characterized by parameters such as speed, angular width, latitude, longitude, and magnetic field flux, among others.

Within the EUHFORIA model, CMEs are inserted into the simulation by assigning the CME’s bulk velocity and plasma properties at the 0.1 au boundary, rather than directly specifying the shock speed V_s. This approach effectively introduces a discontinuity at the inner boundary when the CME enters the computational domain, representing a Riemann problem rather than an initial shock wave (Gogosov 1961; Akhiezer 1975). If the inserted CME’s speed exceeds the local fast magnetosonic speed in the upstream reference frame, a shock wave forms within the computational domain. At the moment of formation, this shock has zero standoff distance from the CME and lacks thus both the shock sheath and the proper lateral extent that would have developed in the solar corona. The shock speed, which will always exceed the CME insertion speed, is determined by the local upstream solar wind conditions and the plasma properties of the CME.

In contrast, to insert a fully resolved fast-mode shock wave with a specific speed (e.g. estimated from white-light images), we start by assuming a shock speed and use the upstream plasma conditions from the EUHFORIA simulation at the inner boundary. The RH ((Rankine 1870; Hugoniot 1887, 1889)) jump conditions for MHD shocks are applied to directly determine the downstream conditions. Assuming a reference frame comoving with the shock, the RH jump conditions can be expressed as

$$\begin{array}{cc}& \left[\rho v_n\right]=0,\hfill \end{array}$$(A.1)

$$\begin{array}{cc}& [(\frac{1}{2}\rho v_n^2+\frac{\gamma }{\gamma -1}p+\frac{B_n^2}{2{\mu }_0})v_n-\frac{\mathit{v}·\mathit{B}}{{\mu }_0}{B}_n]=0,\hfill \end{array}$$(A.2)

$$\begin{array}{cc}& [\rho v_n^2+p+\frac{B_t^2}{2{\mu }_0}]=0,\hfill \end{array}$$(A.3)

$$\begin{array}{cc}& [\rho v_n{\mathit{v}}_t-\frac{B_n}{{\mu }_0}{\mathit{B}}_t]=0,\hfill \end{array}$$(A.4)

$$\begin{array}{cc}& [B_n]=0,\hfill \end{array}$$(A.5)

$$\begin{array}{cc}& [v_n{\mathit{B}}_t-B_n{\mathit{v}}_t]=0,\hfill \end{array}$$(A.6)

where ρ is the plasma density, v is the velocity vector, B is the magnetic field vector, p is the pressure, γ is the adiabatic index, and μ₀ is the permeability of free space. The square brackets [ ⋅ ] denote the jump in the respective quantities across the shock front. The normal vector n is oriented perpendicular to the shock front, and the subscripts ‘n’ and ‘t’ refer to the components of the vectors parallel and perpendicular to the normal direction, respectively.

Solving Eqs. (A.1) to (A.6) provides the downstream plasma conditions, which we prescribe at the inner boundary of the EUHFORIA simulation to model the shock wave with the desired speed. In this work, we focus solely on fast-mode shocks. Additionally, in MHD, the shock compression ratio r = ρ_down/ρ_up approaches (γ + 1)/(γ − 1) for strong shocks, where the fast magnetosonic Mach number M_fms → ∞. For a typical monoatomic gas with γ = 5/3, the compression ratio r tends towards 4. However, EUHFORIA assumes γ = 1.5, resulting in a maximal compression ratio of r → 5. Furthermore, because the magnetic field is almost purely radial at the inner boundary, the shock is what is known as a ‘parallel shock’. In such a configuration and within the assumptions of ideal MHD, the magnetic field has minimal influence on the shock’s behaviour.

To return to the background solar wind state after the shock has passed in the EUHFORIA simulation, we applied the following decay function:

$$\begin{array}{c}\hfill y=y_{s}exp\left[ln\left|\frac{y_b}{y_s}\right|{\left(\frac{t}{\mathrm{\Delta }t}\right)}^{\alpha }\right]\text{with}0\le t\le \mathrm{\Delta }t,\end{array}$$(A.7)

where t denotes the time, y represents any of the plasma variables, y_s denotes the shocked state, and y_b represents the background solar wind state. The parameter Δt defines the time interval over which the plasma variable relaxes to its background state, which we chose in this work as Δt = 20.5 R_⊙/V_s, that is, the time the shock would need to propagate from the solar surface to the inner boundary of EUHFORIA (at 21.5 R_⊙) assuming a constant speed. The parameter α in Eq. (A.7) is a free parameter that controls the decay rate – smaller values of α correspond to faster decays. This decay function is an ad hoc approximation, since a full solution to the downstream profile would require solving the MHD equations from the low corona onwards. Nevertheless, the exponential decay ensures that most of the shocked plasma remains concentrated directly downstream of the shock, qualitatively resembling the SVT blast wave solutions. Additionally, this formulation guarantees a smooth transition from the shock-affected region to the ambient solar wind conditions.

Appendix B: Modified Sedov–von Neumann–Taylor solution

Fig. B.1. Left panel: Upstream number density (n) as a function of shock location (R) along the 0° radial line, with simulation data fitted by n ∝ R−2.05 ± 0.07. Right panel: Downstream pressure immediately behind the shock as a function of R, with simulation data fitted by P ∝ R−3.48 ± 0.18. The symbols indicate in situ spacecraft data.

The SVT solution provides a self-similar description of the propagation of a strong shock wave resulting from an instantaneous, point-like energy release (e.g. explosions) into a homogeneous medium. Due to their simplicity, these solutions have been widely used in astrophysics, notably for describing supernova remnant shock waves. While the interstellar or terrestrial medium where explosions can occur is often assumed to be homogeneous, many astrophysical environments, including the solar wind, exhibit strong density gradients (Cavaliere & Messina 1976; Chevalier 1982).

Once the shock decouples from its CME driver, its evolution becomes similar to a Sedov–Taylor blast wave. At this stage, the shock effectively loses memory of how it originally gained its energy – whether through a sustained piston-driving process or a sudden energy release. Instead, its subsequent evolution is governed by its own momentum and interaction with the ambient solar wind, with its deceleration described well by the Sedov–Taylor framework (e.g. Pinter & Dryer 1990; Moon et al. 2002; Li et al. 2024).

Here, we derive a self-similar solution for a medium where the ambient density decreases with radius as ρ(r) = Ar^−β, with A a proportionality constant. For the expanding solar wind, we can assume β = 2 as a first approximation. To derive the self-similar solution, we assume spherical symmetry, a strong shock with a Mach number ℳ ≫ 1, and energy conservation where the total energy E remains constant over time.

To determine the scaling law for the shock radius R(t) as a function of time t, energy E, and the ambient density parameter A, we used a dimensional analysis similar to that of Sedov (1946). We began by assuming that R(t) has a power-law dependence on time:

$$\begin{array}{c}\hfill R(t)\propto t^{\alpha },\end{array}$$(B.1)

where our objective is to determine the exponent α.

The total energy (E) can be expressed in terms of the post-shock quantities. For dimensional consistency, we have

$$\begin{array}{c}\hfill E\propto \rho (R)R^3{\left(\frac{\mathit{dR}}{\mathit{dt}}\right)}^2.\end{array}$$(B.2)

Substituting $\rho (R)=\frac{A}{R^{\beta }}$ and R ∝ t^α into Eq. (B.2), we obtain

$$\begin{array}{c}\hfill E\propto A{t}^{(5-\beta )\alpha -2}.\end{array}$$(B.3)

Since the total energy E is constant, we solve for α as follows:

$$\begin{array}{c}\hfill \alpha =\frac{2}{5-\beta }.\end{array}$$(B.4)

The shock speed (V(t)) is then given by

$$\begin{array}{c}\hfill V_{\text{sh}}(t)=\frac{\mathit{dR}}{\mathit{dt}}\propto t^{-(3-\beta )/(5-\beta )}.\end{array}$$(B.5)

The pressure immediately downstream of the shock can be estimated using the relation

$$\begin{array}{c}\hfill P\propto \rho V^2,\end{array}$$(B.6)

which follows from the RH momentum jump condition (A.3), assuming negligible upstream magnetic field and upstream pressure (Landau & Lifshitz 1959). Applying the scaling laws for V and ρ, we find

$$\begin{array}{c}\hfill P(t)\propto t^{-6/(5-\beta )}.\end{array}$$(B.7)

In the specific case where β = 2, we have

$$\begin{array}{c}\hfill R(t)\propto t^{2/3},\end{array}$$(B.8)

and thus the shock velocity,

$$\begin{array}{c}\hfill V_{\text{sh}}(t)=\frac{\mathit{dR}}{\mathit{dt}}\propto t^{-1/3}.\end{array}$$(B.9)

Accordingly, the downstream pressure evolves as

$$\begin{array}{c}\hfill P(t)\propto t^{-2}.\end{array}$$(B.10)

Expressing these in terms of the shock radius R, using t ∝ R^1/α, we obtain

$$\begin{array}{c}\hfill V_{\text{sh}}(R)\propto R^{-1/2},\end{array}$$(B.11)

and

$$\begin{array}{c}\hfill P(R)\propto R^{-3}.\end{array}$$(B.12)

In summary, we find that the shock radius expands as t^2/3, indicating a decelerating expansion. This leads to a reduction in shock velocity over time as t^−1/3 or with distance as R^−1/2, consistent with, for example, Cavaliere & Messina (1976) and Smart & Shea (1985). Meanwhile, the downstream pressure diminishes rapidly, scaling as t⁻² or R⁻³.

These theoretical predictions align well with the observed trends shown in Fig. B.1, where the fitted power laws are n ∝ R^{−2.05 ± 0.07} for the upstream density and P ∝ R^{−3.58 ± 0.18} for the downstream pressure. While the power laws capture the average decay rates approximately, significant deviations are introduced by the preceding CMEs, which create fluctuations around these mean trends.

In addition, in situ measurements from various spacecraft are included in Fig. B.1. It is important to note, however, that none of the spacecraft are exactly positioned along the 0° line in the solar equatorial plane. Earth, for instance, is closest to this line but is situated approximately 7 degrees south of the equatorial plane. Despite these positional discrepancies, the in situ data generally exhibit a similar decreasing trend, particularly in density.

For the downstream thermal pressure, both the simulation and observations display a comparable decaying trend; however, the simulations consistently yield higher pressure values than those observed. This difference is also evident in Fig. 3, which compares the simulated and observed pressure time profiles. The likely causes for this discrepancy include a slight overestimation of the shock strength in the simulations, as well as an underestimation of the magnetic field magnitude, leading to a lower magnetic pressure in the simulation compared to actual observations.

All Tables

All Figures

Fig. 1. Spacecraft constellation in the solar equatorial plane, showing nominal Parker spiral field lines connecting each spacecraft to the Sun. The field lines were derived using measured solar wind speeds for STEREO-A, Earth, and Parker Solar Probe, while for Solar Orbiter and BepiColombo, where no solar wind plasma data were available at the time of the event, solar wind speeds from the EUHFORIA simulation are used. The two arrows indicate the estimated locations of the ARs from which the CMEs erupted. The dashed pink line represents the combined extent of the two shock waves injected in the EUHFORIA simulation. This figure was made using the Solar-Mach tool (Gieseler et al. 2023).

In the text

Fig. 2. Selected snapshots of the CME evolution on 13 March 2023, as seen in coronagraphic observations from LASCO and STEREO-A. Each panel shows an overlay of a plain coronagraphic image and the corresponding running difference image at the same time from LASCO-C2 (panels a and b) and COR2-A (panels c to f). At selected frames we note the CMEs and outline the location of the CME-driven shock wave. The interaction of the shock with streamers (S1 and S2) are also noted in panels d to f (streamers deflections S1D and S2D). The various features are discussed in detail in Sect. 2.

In the text

Fig. 3. Comparison of observed (blue) and simulated (orange) solar wind parameters at various spacecraft locations: Parker Solar Probe (top left), Solar Orbiter (top right), L1 (ACE and WIND; bottom left), and STEREO-A (bottom right). For each spacecraft, the panels display: (a) solar wind speed, (b) proton number density, (c) thermal proton pressure, (d) magnetic field magnitude, (e) magnetic latitudinal angle, and (f) magnetic azimuthal angle.

In the text

	Fig. 4. Snapshot of the EUHFORIA simulation on 14 March at 05:57 UT. Left panels: Solar equatorial plane. Right panels: Meridional slices containing Earth. Top row: Solar wind speed. Bottom row: Divergence of the solar wind velocity ∇ ⋅ Vsw. Dashed lines indicate projections of magnetic field lines connecting to various spacecraft. A video version of this figure is available online.
In the text

Fig. 5. Main-event shock surface, along with the magnetic field lines connecting Earth (green), STEREO-A (red), Solar Orbiter (blue), BepiColombo (orange), and Parker Solar Probe (purple). In the top panel, the shock surface is colour-coded by the fast magnetosonic Mach number (Mfms), while in the bottom panel, it is colour-coded by the shock obliquity angle (θBn). The associated online video provides a 360° view around the shock surface.

In the text

Fig. 6. Evolution of the main-event shock flanks along two radial directions in the solar equatorial plane. The blue symbols represent the line at 0°, aligned with the positive x-axis in Fig. 4, while the green symbols represent a line rotated 135° anti-clockwise from the positive x-axis. The panels present: (a) hours elapsed since the blast wave injection, including power-law fits and observed spacecraft shock crossing times; (b) shock speed; (c) upstream solar wind speed; (d) upstream solar wind density, scaled with r2; (e) upstream fast-magnetosonic speed; and (f) fast-magnetosonic Mach number (Mfms).

In the text

	Fig. B.1. Left panel: Upstream number density (n) as a function of shock location (R) along the 0° radial line, with simulation data fitted by n ∝ R−2.05 ± 0.07. Right panel: Downstream pressure immediately behind the shock as a function of R, with simulation data fitted by P ∝ R−3.48 ± 0.18. The symbols indicate in situ spacecraft data.
In the text