Hybrid-Vlasov modeling of three-dimensional dayside magnetopause reconnection

Dayside magnetic reconnection at the magnetopause, which is a major driver of space weather, is studied for the first time in a threedimensional (3D) realistic setup using a hybrid-Vlasov kinetic model. A noon–midnight meridional plane simulation is extended in the dawn–dusk direction to cover 7 Earth radii. The southward interplanetary magnetic field causes magnetic reconnection to occur at the subsolar magnetopause. Perturbations arising from kinetic instabilities in the magnetosheath appear to modulate the reconnection. Its characteristics are consistent with multiple, bursty, and patchy magnetopause reconnection. It is shown that the kinetic behavior of the plasma, as simulated by the model, has consequences on the applicability of methods such as the four-field junction to identify and analyze magnetic reconnection in 3D kinetic simulations. VC 2020 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (http:// creativecommons.org/licenses/by/4.0/). https://doi.org/10.1063/5.0020685


I. INTRODUCTION
Dayside magnetic reconnection is a major driver of magnetospheric dynamics (Dungey, 1961), allowing energy from the solar wind to penetrate the magnetosphere and to cause some of the most spectacular effects of space weather like aurora. It can occur at the Earth's dayside magnetopause under a variety of conditions, which affect the location, shape, and dynamics of the reconnection sites.
The two main modes of dayside magnetic reconnection that have been proposed and observed are antiparallel reconnection, where the two reconnecting magnetic fields form an angle close to 180 (e.g., Crooker, 1979;Luhmann et al., 1984), and component reconnection, where oppositely directed components reconnect in the presence of a significant guide field (e.g., Sonnerup, 1974;Gonzalez and Mozer, 1974). They are combined in the maximum magnetic shear model, giving an estimate of the dayside magnetic reconnection line location depending on interplanetary magnetic field (IMF) conditions (Trattner et al., 2007;Trattner et al., 2012Trattner et al., , 2017Trattner et al., , 2018. This model, built on in situ and remote space-and groundbased observations, successfully predicts the position of dayside magnetic reconnection in a statistical sense. It does not, however, describe all the local and dynamical properties that have been observed (see, e.g., the review by Fuselier and Lewis, 2011, and references therein).
In support of inherently sparse spacecraft and ground-based observations, there is a long history of studying dayside reconnection in global simulations using magnetohydrodynamic (MHD) models in 2D and 3D setups (e.g., Ogino et al., 1986;Shi, Wu, and Lee, 1988;Laitinen et al., 2006;Borovsky et al., 2008;Hoilijoki et al., 2014;Komar et al., 2015). Efforts have been undertaken using kinetic models as well, although their computational cost remains an issue (Hesse and Cassak, 2020). Therefore, approaches such as propagating particles or embedding kinetic regions in global MHD models, reducing the simulations to two-dimensional (2D) slices, and/or scaling the magnetic dipole strength are required (e.g., Nakamura and Scholer, 2000;Karimabadi et al., 2014;Connor et al., 2015;Peng et al., 2015;Chen et al., 2017;T oth et al., 2017;Walker et al., 2019). There is still, however, no global, kinetic, self-consistent, and realistically scaled model describing the high spatial and temporal variability of dayside magnetopause reconnection that has been abundantly observed.
The hybrid-Vlasov model, Vlasiator , simulates kinetic ions coupled to fluid electrons in near-Earth space. It has been used to study dayside magnetopause reconnection, but only in 2D noon-midnight meridional planar simulations so far, due to the prohibitive cost of extending the global domain into 3D. These 2D investigations have uncovered large-scale perturbations of the magnetosheath and the bow shock related to reconnection Jarvinen et al., 2018), described the impact dayside reconnection can have on tail reconnection , confirmed analytic predictions of asymmetric dayside reconnection rates , and described the dynamics of reconnection sites and the properties of the FTEs generated at the dayside magnetopause (Hoilijoki et al., , 2019Akhavan-Tafti et al., 2020). The present study shows the results from a Vlasiator simulation where, for the first time, the system is extended in the third dimension to cover 7 R E in the dawn-dusk direction. This causes the magnetosphere-like structure to assume a cylindrical geometry, allowing to retain periodic boundaries in the dawn-dusk direction. The resulting numerical experiment confirms the findings of Hoilijoki et al. (2017) that mirror mode waves present in the magnetosheath affect dayside reconnection also in 3D. This novel setup further allows us to investigate the high temporal and spatial variability of dayside magnetic reconnection, occurring in the presence of a steady southward IMF.
Section II presents the simulation and analysis methods used, the results are presented in Sec. III and discussed in Sec. IV, and a summary and conclusions are given in Sec. V.

II. METHODS A. Hybrid-Vlasov simulation
Vlasiator is a hybrid-Vlasov simulation geared toward modeling the coupled solar wind-magnetosphere-ionosphere system in 3D, using the ion-kinetic hybrid approximation. It solves Vlasov's equation for ions coupled to Maxwell's equations closed with a generalized Ohm law, including the Hall term. This results in a kinetic description of the proton population while approximating the electrons as a cold and massless charge-neutralizing fluid (von Alfthan et al., 2014;Palmroth et al., 2018;Vlasiator Team, 2020).
All previous self-consistent magnetospheric Vlasiator studies used 2D setups. For simulations in the noon-midnight meridional plane, a line dipole is used with a scaling so as to yield a realistic magnetopause standoff distance (Daldorff et al., 2014). In the present work, the simulation plane is extended in the dawn-dusk direction, yielding a 3D simulation box. The line dipole is still used so that the magnetosphere-like structure assumes a cylindrical geometry under periodic boundary conditions in the dawn-dusk direction. This numerical experiment allows direct comparison to previous work in 2D (Hoilijoki et al., , 2019Akhavan-Tafti et al., 2020) at a fraction of the computational cost of a full 3D simulation, including the 3D geomagnetic point dipole.
The upstream fast solar wind conditions are steady with a Maxwellian velocity distribution of density 1 cm -3 , temperature 0.5 MK, and velocity -750 km s À1 along X, as well as a purely southward IMF of B Z ¼ À5 nT strength. Hence, the solar wind ion inertial length is 227 km and the ion thermal gyroradius is 134 km. The simulation box spans X 2 ðÀ16; 31Þ R E ; Y 2 ðÀ3:5; 3:5Þ R E ; Z 2 ðÀ35; 35Þ R E in the Geocentric Solar-Ecliptic (GSE) coordinate system (without a dipole tilt), with 200, 30, and 300 cells, respectively, yielding a spatial resolution of 0.24 R E . This spatial resolution is approximately 2-5 times coarser than that used in 2D studies with Vlasiator, but it is comparable to the finest resolutions used typically in MHD modeling. For example, T oth et al. 2 R E in their high-resolution study. In comparison, the resulting ion skin depth is of the order of 100 km or 0.016 R E in the magnetosheath and at the magnetopause. The velocity resolution is 50.25 km/s, which is sufficient, given the spatial resolution . The phase space density sparsity threshold, below which phase space cells are disregarded and not propagated, is applied in this run . It is set to 10 À17 s 3 m -6 , where the plasma density is below 0.01 cm -3 and 10 À15 s 3 m -6 above 0.1 cm -3 , with a linear ramp between those values. The -X and 6Z walls have a Neumann condition ensuring outflow of the plasma, whereas the Y direction is treated periodically. The perfectly conducting inner boundary is cylindrical with a radius of 4.7 R E and has a constant and static Maxwellian velocity distribution of the same density and temperature as the solar wind. This choice of distribution has proven to typically provide better stability during the initialization of Vlasiator runs. The initial setup consists of the line dipole field as well as the solar wind plasma permeating the whole box. The southward IMF flows into the box from the sunward þX boundary. The simulation runs for a total of 761.5 s and the southward IMF reaches the dayside magnetopause before t ¼ 650 s. The run took around 0.3 Â 10 6 core-hours to perform and the whole dataset takes 1 TB on disk.

B. Identification of dayside reconnection
During the southward IMF, dayside reconnection occurs in the subsolar region. It can be searched for using the four-field junction method, which was developed for global MHD simulations, and relies on the magnetic field topology. In a volume covering the dayside magnetosphere, magnetic field lines are traced forward and backward starting from every simulation cell. The magnetic connection of each cell to the outer walls (open; solar wind) or the inner boundary (closed; ionosphere) is recorded. There thus exist four types of field lines, namely those open at both ends, open-closed, closed-open, and closed at both ends. Magnetic reconnection, which connects volumes threaded by magnetic fields which were previously distinct, can occur in the vicinity of any point near which all four types of magnetic field lines are present simultaneously (Laitinen et al., 2006). In this study, the nearest (3 Â 3 Â 3 cells) and second-nearest (5 Â 5 Â 5 cells) neighbors to each cell are taken into account to determine the location of the four-field junction. Active magnetic reconnection near fourfield junction points can then be identified based on the ongoing topological reconfiguration of magnetic fields and energy conversion, evident for example in the strong acceleration of plasma in the exhaust regions.

ARTICLE
scitation.org/journal/php It is, however, possible that magnetic reconnection occurs between domains not comprising the four types of field lines mentioned above. A typical example are FTEs with active reconnection regions on either side of a flux rope. The flux rope magnetic field lines do not necessarily fulfill the four-field junction requirements for both sides simultaneously. Another example is lobe reconnection beyond the cusps during the northward IMF. In such cases, careful analysis is required in order to identify reconnection sites based on the magnetic field configuration and the signatures of the energy conversion inherent to magnetic reconnection.  Figure 1 shows the magnitude of the magnetic field B and Fig. 2 shows the proton number density n. The bow shock forms the interface between the solar wind flowing in from the þX side and the compressed, heated plasma in the magnetosheath. The magnetic field magnitude B (Fig. 1) shows well the magnetopause in panels (a)-(d), as the dominant component B Z goes to zero at the transition between the magnetosheath, where B is southward, and the magnetosphere, where the dipolar B is northward. Both Figs. 1 and 2 also include a pair of contours of the north-south component of velocity at V Z ¼ 6300 km/s. These contours delineate narrow high-velocity channels of plasma escaping along the magnetopause from the dayside toward and beyond the polar cusp region. The channels are the exhausts of active magnetic reconnection between the southward magnetosheath field and the northward dipole field. They are locally super-Alfv enic and the velocity corresponds to the expected order of magnitude considering the inflow densities and magnetic fields (Cassak and Shay, 2007).

A. Magnetosheath structure
The upstream solar wind is constant and uniform by design in this simulation and the line dipole magnetic field has a translation symmetry in the Y direction. Furthermore, neither the line dipole field nor the IMF have a B Y component. However, the cuts at a constant X in Figs. 1 and 2 show that some structures in the magnetosheath break the translation symmetry along Y. At X ¼ 16 R E [panels (f)], the perturbations are limited albeit already present between Z ¼ 0 and 5 R E . At X ¼ 12:5 R E [panels (e)], perturbations along Y are more widespread, especially in the southern half. In the planes cutting the magnetopause [panels (b)-(d)], the magnetopause transition is highlighted by the drop in B and partially by the V Z contours. The magnetopause is clearly not straight and parallel to the Y direction, unlike what would be expected from the interaction of the IMF and the line dipole in the absence of perturbations. A range of wave mode numbers m can be identified in the perturbations in the Y direction. In Figs. 1(c) and 2(c), the blue contour near Z ¼ À5 R E has a wave mode number m ¼ 1. One full sine-like oscillation m ¼ 2 is visible, for example, in Figs. 1(f) and 2(f) near Z ¼ 3 R E and in panels (b) and (c) near Z ¼ À7 R E , m ¼ 4 (two full oscillations) is visible, for example, near the northern magnetopause intersection in panels (b-d) and Figs. 2

B. Magnetopause reconnection: Four-field junction
The first step to study the location of dayside magnetic reconnection in this setup is to use the four-field junction method described in Sec. II, which has a heritage in the analysis of global MHD simulations of the magnetosphere and forms the basis of more efficient identification algorithms (e.g., Komar et al., 2013;Glocer et al., 2016). Figure 3 shows the isosurface B Z ¼ 0, at which the southward magnetosheath magnetic field flips to the northward line dipole field, corresponding to the dayside magnetopause in the present setup. The color code shows the V Z component of the plasma velocity, as the plasma flow is mainly expected to diverge symmetrically to the north and south of the Z ¼ 0 plane, accelerated by magnetic reconnection. While the flow diverges globally, as constrained by the geometric setup, the structure shown in the non-saturated region of the colorcoded isosurface is complex. Several spots of southward flow embedded in northward flow and vice versa are present in the region where jZj < 1 R E . The magnetic field reversal and the center of the reconnection exhaust jets are not necessarily colocated in asymmetric reconnection (e.g., Cassak and Shay, 2007), so that the displayed V Z is not the   . The reconnection line, as represented by the set of identified four-field junction points, lies between À1 R E > Z > À2:5 R E in Fig. 3. It can be very narrow, less than 0.5 R E , and near Y ¼ 0, but it thickens to 1.5 R E near the periodic edge of the simulation box. Over the time covered by Animation 1 (Multimedia view), the thickness of the line along Z reaches up to 3.1 R E , while its maximum southward (northward) extent is -3.6 R E (0.7 R E ). It can also become disjoint, like at times around 698 or 714 s. Figure 4 shows the simulation at the same time as Fig. 3, but from a rotated perspective. The isosurface corresponds again to the magnetopause where B Z ¼ 0 and the color corresponds to the V Z component. A selection of magnetic field lines highlights several reconnection sites, which were identified by visual inspection. As expected, the four-field junction (shown only for the points upstream of the B Z ¼ 0 isosurface) coincides with an X line configuration, as demonstrated in the south near Y ¼ 1:5 R E . Three other reconnection regions are shown further north near Y ¼ À2:2; 0; and 1 R E . They are not associated with the junction of the four field line types, as the southern field line examples are likely not connected to the closed boundary due to the reconnection site further south. They however coincide with regions of north-south flow divergence.

C. Beyond the four-field junction
The sign of the flow component V Z is compared between the north and south neighbors of every simulation cell straddling the magnetopause surface. If V Z has a different sign on either side, the given cell is identified as a candidate location for magnetic reconnection. The time series of the Z coordinate for all points identified with the four-field junction and the V Z divergence methods is shown in Fig.  5(a). Figure 5(b) additionally shows the histograms of the three sets of points. The V Z divergence method is only applied at the magnetopause, whereas four-field junction points can be identified at multiple X coordinates for a given pair of (Y, Z) coordinates. Therefore, the four-field junction histograms are built from the sets of unique (Y, Z) coordinates, disregarding the X coordinate. Figure 5 shows that there are broad regions, where either method detects points individually as well as regions where both methods detect points simultaneously. The histograms are skewed southward, a feature discussed in more detail below.

IV. DISCUSSION A. Magnetosheath waves and symmetry breaking
Previous detailed analysis of 2D Vlasiator simulations shows that mirror mode waves feature prominently in the magnetosheath . The comparison of Figs. 1 and 2, helped by their similar color ranges, demonstrates that all the large-scale magnetosheath wave-like structures in B and n are anticorrelated, in the subsolar magnetosheath [panels (a), (e), and (f)] as well as further down the flanks [panel (a)]. This is characteristic of mirror mode waves, which are thus dominating the magnetosheath in the present simulation.  Fig. 3 in a perspective view, including a selection of magnetic field lines highlighting magnetic reconnection sites. The four-field junction points are only shown upstream of the B Z ¼ 0 isosurface. The field lines are colored by the X coordinate and dotted where they are behind the isosurface. Note the stretching of the X axis to ease the view, which exaggerates the apparent curvature of the isosurface.

ARTICLE
scitation.org/journal/php Dubart et al. (2020) investigated the impact spatial resolution has on the growth of magnetosheath waves in 2D Vlasiator simulations. They show that, while at low resolutions the electromagnetic ioncyclotron (EMIC) waves are almost absent, the mirror mode is still present. Faint stripes near the bow shock at X < 7 R E and jZj > 23 R E show a correlation between B and n, indicative of EMIC waves. They are clearly underresolved at such wavelengths close to the spatial grid resolution, meaning also that the mirror modes are over-emphasized, as less energy is dissipated by the EMIC waves.
In this 3D setup, the dawn-dusk Y direction covers 7 R E and 30 grid points, and the upstream and inner boundary conditions enforce constant parameters along the Y direction. The results presented in Sec. III A demonstrate that in spite of these constraints, the instabilities growing in the magnetosheath break the translation symmetry in the Y direction and non-uniformities with mode numbers between m ¼ 1 and m ¼ 6 build up. No other source of perturbations seems to be present, therefore this magnetosheath activity is most likely what breaks the symmetry and in the end affects the shape and dynamics of the magnetopause and of the magnetic reconnection described in Sec. III.

B. Occurrence of magnetic reconnection
Vlasiator does not include an explicit resistive term Vlasiator Team, 2020), hence magnetic reconnection is made possible by the numerical diffusivity or resistivity arising from the finite order of accuracy of the algorithms, as is the case in global simulations without such an explicit term Palmroth et al., 2018). The Hall term in Ohm's law ensures a fast reconnection rate (Birn et al., 2001). Numerical resistivity is expected to depend on the grid resolution and maybe on the orientation of the grid with respect to the current sheet. However, a convergence study or the comparison to a run with an explicit resistive term is outside the scope of the present work. It should be noted, however, that the numerical diffusivity is low enough so as to not adversely impact the modeled structures. This can be assessed in particular in Figs. 1 and 2, where it is evident that the bow shock jump and the magnetopause magnetic field reversal are resolved sharply at the grid resolution.
Dayside reconnection in 2D Vlasiator simulations has been shown to proceed at rates consistent with the analytic model of Cassak and Shay (2007)  . Whereas in 2D it is straightforward to map the reconnection rate to the out-of-plane component of the electric field, determining the direction of that component is not as easy in 3D, given the shape the reconnection line can assume. Going back to the formulas of Cassak and Shay (2007), picking upstream values for B and n can be done in a direction orthogonal to the magnetopause. However, it remains difficult to determine systematically the outflow direction and the distance at which the outflow values should be taken to compute an estimate of the reconnection rate at every four-field junction point. The flow pattern in the region where jZj < 2 R E in Figs. 3 and 4, resulting from the interaction of magnetosheath perturbations and reconnection exhausts oriented in multiple directions, makes the systematic determination of an outflow direction nearly impossible. Therefore, an estimation of the local reconnection rate has not been undertaken in the present setup.

C. Identification of magnetic reconnection
The variability of the location, velocity, and reconnection rate of the X points forming at the dayside magnetopause under a southward IMF has been demonstrated in previous work based on 2D simulations with Vlasiator (Hoilijoki et al., , 2019Akhavan-Tafti et al., 2020). There are convenient methods used to identify X and O points (reconnection sites and magnetic islands) in 2D, based on the local behavior of the magnetic flux function (Servidio et al., 2010). They are, however, not applicable in 3D without significant constraints incompatible with the present setup (Yeates and Hornig, 2011).
The four-field junction is inherently built upon the global topological properties of the magnetosphere, in contrast to the local

ARTICLE
scitation.org/journal/php methods based on the magnetic flux function in 2D. More general methods than the four-field junction exist in 3D, which do not rely on the connectivity of the field lines with the ionosphere. In particular, the identification of the magnetic skeleton of a given magnetic field configuration, comprising magnetic null points, spines, separatrix surfaces, and separators, is a powerful tool since magnetic reconnection takes place along separators joining pairs of null points (Haynes and Parnell, 2010). The first step in determining the magnetic skeleton is finding all the magnetic null points, for which methods have been devised both for simulations (e.g., Haynes and Parnell, 2007) and observations (e.g., Fu et al., 2015). In general though, algorithms such as the trilinear method proposed by Haynes and Parnell (2007) call for sufficient spatial resolution of the magnetic field and its derivatives. In the present study, the resolution is not fine enough to yield reliable interpolation results in the vicinity of reconnection points. Moreover, a systematic study of magnetic nulls in a range of space plasma configurations simulated with the iPIC3D code shows another limit affecting kinetic simulations (Olshevsky et al., 2016). The authors identify more than 8000 magnetic nulls in their dipolar mini-magnetosphere setup (Figure and Section 6, Olshevsky et al., 2016). This makes a systematic determination of separators between pairs of nulls like in MHD simulations (e.g., Haynes and Parnell, 2010;Glocer et al., 2016) impractical or simply impossible in kinetic simulations. The present work shows that the global four-field junction method does not detect magnetic reconnection sites comprehensively, as kinetic effects causing inhomogeneous perturbations trigger multiple reconnection sites as shown in Fig. 4. Topologically, a field line passing by two distinct reconnection sites cannot be registered through four-field identification at both sites. Multiple reconnection is of course a generation mechanism of flux ropes and FTEs, one of which is located in the region partially enclosed by the two rightmost reconnection sites in Fig. 4 and another one southward of the leftmost reconnection site. Consequently, the four-field junction method is insufficient for analyzing the structure and dynamics of magnetic reconnection sites systematically and automatically in a 3D kinetic simulation, such as the one presented here. Figure 5 compares the four-field junction location with the sites, where V Z flips directions. First, it should be noted that the divergence of the flow is dictated by fluid dynamics even in the absence of magnetic reconnection. Nevertheless, the presence of strong, diverging jets is an indicator of possible active reconnection. Second, it should also be noted that the reversal of V Z is not expected to be seen at higher jZj, as the reconnection sites are likely to be advecting along with the magnetosheath flow accelerating toward the flanks. Indeed, it is not visible beyond Z > 1:5 R E (Z < À2:5 R E ) during the time covered by Fig. 5 and Animation 1 (Multimedia view). The resolution in this study does not allow to evaluate reliably the derivatives of the flow at the magnetopause with respect to the ambient magnetosheath flow, in order to track diverging flows at higher jZj. While there is some overlap between regions where the four-field junction detects reconnection and regions of V Z reversal, it is partial and broad areas covered by only one of the regions exist as well. This indicates that the analysis of the V Z direction is not enough to complement the four-field junction method in detecting magnetic reconnection. Furthermore, it has been shown that the flow stagnation point and the magnetic reconnection point (X point in 2D) are not colocated in asymmetric reconnection (e.g., Cassak and Shay, 2007), but the difference in position is not expected to be large compared to the resolution discussed in this study.
Animation 1 (Multimedia view) and Fig. 5 demonstrate a bias locating the four-field junction points mostly south of the equatorial plane, despite the symmetry of the boundary conditions. However, the V Z flip locations shown in Fig. 5 and the example reconnection sites highlighted by magnetic field lines shown in Fig. 4 are more evenly distributed around Z ¼ 0. This illustrates yet further that the four-field junction method does not detect all reconnection sites, leading to potentially erroneous conclusions in the absence of further analyses.
All of this shows that both the topology of the magnetic field and the detailed dynamics of energy conversion have to be considered locally in order to identify magnetic reconnection and its effects. A robust, automatized algorithm detecting sites of magnetic reconnection in 3D global kinetic simulations would be a valuable tool and remains the topic of current research.

D. Characterization of dayside magnetic reconnection
The first statistical and geometrical expectation might be to find the reconnection line at Z ¼ 0 [see, e.g., Fig. 4.2(c) by Hoilijoki et al., 2014, based on MHD modeling]. But in a configuration where the IMF is purely southward and not tilted with respect to the dipole axis, the region where the dipole magnetic field is nearly antiparallel to the magnetosheath field covers a large fraction of the subsolar magnetopause (see Fig. 7 by Trattner et al., 2012, for example). The antiparallel region would even cover the whole magnetopause sunward of the cusps in the absence of inhomogeneities in the 3D setup presented here. In the presence of the said inhomogeneities, it is therefore not surprising to find the locus of dayside magnetic reconnection fluctuating by several R E in Z, showing a complex structure and sometimes even a disjoint topology.
Although the simulation presented here covers 7 R E in Y and lacks the dawn-dusk curvature of the 3D geomagnetic dipole, its results are consistent with several decades of observations of multiple, bursty, and patchy dayside magnetopause reconnection. Multiple reconnection occurring simultaneously at different Z is obvious in Fig. 4, while Animation 1 (Multimedia view) shows that reconnection sites can be activated at multiple Z and then travel in Z. Multiple reconnection as well as bursty (or intermittent) reconnection have long been invoked as the generation mechanism for FTEs (e.g., Russell and Elphic, 1978;Kan, 1988;Fear et al., 2008Fear et al., , 2009Eastwood et al., 2016). Patchy reconnection is prevalent in the present simulation: throughout the time covered by Animation 1 (Multimedia view), the four-field junction points can be divided into distinct groups, constrained to sizes of a few R E by the simulation dimensions. Patchy reconnection has been proposed and investigated theoretically (e.g., Galeev, Kuznetsova, and Zeleny, 1986;Kan, 1988;La Belle-Hamer, Fu, and Lee, 1988;Kan et al., 1996), and subsequently observed in situ and remotely at the dayside magnetopause (e.g., Milan et al., 2016;Walsh et al., 2017;Zou et al., 2019Zou et al., , 2020. Of course, the present simulation results do not preclude the existence of extended, continuous reconnection lines, which have also been observed (e.g., Phan et al., 2004;Milan et al., 2016;Zou et al., 2019) and would likely also occur in larger, fully 3D simulations.

V. CONCLUSIONS
We use a hybrid-Vlasov simulation to model dayside magnetopause reconnection in a three-dimensional setup by extending the typical noon-midnight meridional plane domain in the dawn-dusk direction. Despite the constant southward interplanetary magnetic field, the uniform inner boundary condition, and the relatively coarse spatial resolution, kinetic instabilities in the magnetosheath break the expected symmetries so that the site of dayside reconnection travels north and south of the equatorial plane and it is neither a straight nor a continuous line. The results are consistent with multiple reconnection occurring simultaneously at several latitudes, with bursty or intermittent reconnection varying in time, as well as with patchy reconnection localized at small spatial scales. Both the four-field junction method and the flip of the north-south flow component are used to identify some locations of magnetic reconnection. We show, however, that these methods need to be complemented as further reconnection sites not matching either of them were identified in the simulation run. Systematic methods achieving the goal of detecting all magnetic reconnection sites in global kinetic simulations remain under investigation and will be the topic of further research.

ACKNOWLEDGMENTS
We acknowledge the European Research Council for Starting Grant No. 200141-QuESpace, with which Vlasiator Vlasiator Team, 2020)  Computing resources provided by the CSC-IT Center for Science in Finland were used to produce the results presented.

DATA AVAILABILITY
The data of the run are accessible by following the rules published at http://www.physics.helsinki.fi/vlasiator as a consequence of the large dataset size.