Effects of dispersion of the dust velocity in the LISM on the interstellar dust distribution inside the heliosphere
EEffects of dispersion of the dust velocity in the LISM on theinterstellar dust distribution inside the heliosphere
E.A.Godenko , , ∗ , V.V.Izmodenov , , Institute for Problems in Mechanics RAS, 119526, Moscow, Pr.Vernadskogo, 101-1 Lomonosov MSU, Moscow Center for Fundamental and Applied Mathematics, 119992,Moscow, GSP-1 Leninskie Gory Space Research Institute RAS, 117997, Moscow, Profsoyuznaya Str 84/32 Annotation – Interstellar dust (ISD) penetrates into the heliosphere due to the relativemotion of the Sun and the local interstellar medium (LISM). Inside the heliosphere and atthe boundaries, where solar wind interacts with the LISM, distribution of ISD is modifieddue to the action of the electromagnetic forces, the solar gravitation and the radiationpressure. These forces make the distribution of the ISD particles in the heliosphere inho-mogeneous. In previous work we demonstrated the existence of singularities in the ISDdensity distribution at 0.03 - 10 a.u. north and south with respect to the heliospheric cur-rent sheet. In this paper we show that dispersion in the ISD velocity distribution stronglyaffects the singularities. Even small values of dispersion have the drastic impact on thedensity distribution and smooth the high density layers discovered previously.
Key words: dust, heliosphere, numerical methods. ∗ E-mail < [email protected] > a r X i v : . [ a s t r o - ph . S R ] F e b INTRODUCTION
The local interstellar medium (LISM) moves relative to the Sun with the speed ∼
26 km/s(Witte 2004, McComas 2015). Besides the plasma and neutral components, the LISM alsocontains dust component (Mann 2010). Unlike the plasma particles, the neutral and dustparticles can penetrate into the heliosphere due to the relative motion. For example, themean free path of neutral hydrogen due to charge exchange is ∼ −
100 a.u. (Izmodenovet al. 2000), comparable with the characteristic size of the heliosphere.The first evidence for the existence of interstellar particles in the heliosphere was theLyman- α emission of interstellar neutrals (Bertaux, Blamont 1971). Direct measurementsof the interstellar helium atoms were obtained by the Ulysses/GAS instrument (Witte1992), and since the mean free path of interstellar helium is much larger than the sizeof the heliosphere, one can derive the macroscopic parameters of the LISM from thesemeasurements. Nowadays, direct measurements of interstellar neutrals (hydrogen, oxygenand helium) are performed on IBEX using the IBEX-Lo instrument (e.g. Moebius etal. 2009, Katushkina et al. 2015, Baliukin et al. 2017). On the spacecraft SOHO(SWAN instrument) measurements of intensity and spectral characteristics of the Lyman- α emission are continuing (e.g. Qu´emerais et al. 2013). Various models of the heliosphereare employed for the analysis of the experimental data (e.g. Izmodenov, Alexashov 2015,2020, Pogorelov et al. 2011, Zirnstein et al. 2016).The ISD grains are solid grains with characteristic sizes in the range of hundreds of nanome-ters to microns (Mathis et al. 1977). Chemical composition of ISD is carbonaceous ma-terials and astronomical silicates (Draine 2009). The mass fraction of ISD in the LISM isabout 1 % (Mann 2010). The ISD grains are charged positively as net effect of differentphysical processes such as photoelectron and secondary electron emissions. The presenceof nonzero electric charge makes the trajectories of ISD more complex than of interstellarneutrals (not taking charge exchange with protons into account).It is difficult to detect ISD in the heliopshere because of the presence of interplanetary dust,which is emitted from asteroids, comets and other large objects in the Solar system. It isgenerally supposed that in the undisturbed LISM the interstellar dust is comoving withother components. This assumption was used in order to detect the ISD grains on Ulysses(Gr¨un et al. 1994). Moreover, the trajectory of Ulysses went significantly out of the eclipticplane and thus it gave an opportunity to relatively easily separate interstellar dust frominterplanetary dust, which is located principally in the ecliptic plane (e.g. zodiacal dust).Presence of ISD was also confirmed in the measurements on board the Galileo (Altobelliet al. 2005) and Cassini (Altobelli et al. 2007) spacecraft.The first models of the ISD distribution in the heliosphere were made by Bertaux, Blamont(1976) and Levy, Jokipii (1976). They studied the distinct influence of the gravitationaland electromagnetic forces on the motion of the dust particles in the heliosphere. The nextwave of interest in the ISD studying was associated with Ulysses measurements. Landgrafet al. (2000, 2003) analyzed these measurements using the Monte-Carlo modeling. Theyconsidered the combined influence of the gravitational, radiation pressure and electromag- 3 –netic forces on the particles in presence of time-dependent solar magnetic field. The ISDdistribution and filtration of dust grains by the magnetic field at the heliospheric bound-aries were explored by Czechowski, Mann (2003), Alexashov et al. (2016). Slavin et al.(2012) have built a 3D model of the ISD distribution for two opposite phases of the helio-spheric magnetic field (focusing and defocusing). It is also taken account of the turbulenceof the interstellar magnetic field and dependence of the surface charge potential on theheliocentric distance. Nowadays, the Monte-Carlo method is often used for theoreticalstudies of ISD. The descriptions and results of the modeling are shown in Sterken at al.(2012, 2019), Strub et al. (2015, 2019). These models are developed from the earlier modelof Landgraf et al. (2000) using advanced numerical techniques and taking into accountof the newer measurements. Mishchenko et al. (2020) applied a Lagrangian method (seeOsiptsov 2000) to discover singularities in the distribution of ISD in the heliosphere. In thesimplified stationary case when the heliospheric current sheet is a plane coinciding withthe solar equatorial plane they demonstrated the existence of density singularities wherethe number density is infinite. They showed that the singularities form several dense dustlayers for each size of the ISD particles on both sides of the current sheet. These singu-larities have never been observed in the previous papers studying dust distribution in theheliosphere by Monte-Carlo modeling because it requires a computational grid with anextremely high spatial resolution. In this paper we use a computational grid with cell sizeof 10 − a.u. and for studying local effects near density peculiarities — of 10 − a.u.Mishchenko et al. (2020) used the assumption that the ISD particles have identical veloci-ties in the LISM. Due to the fact that the ISD particles have nonzero electric charge, theyinteract with the interstellar magnetic field. Fluctuations of the magnetic field lead to theacceleration of the charged dust particles (Hoang et al. 2012), that breaks the uniformityin the ISD velocity distribution in the LISM and adds some rather small dispersion. Thegoal of this paper is to study the influence of dispersion in the velocity distribution ofISD in the LISM on the emergence of the density singularities in the heliosphere. Slavinet al. (2012) also explore dispersion in the undisturbed LISM, but they do not study itsinfluence on the singularities, since the computational grid used is quite coarse (5 a.u. foreach direction). DESCRIPTION OF THE MODEL
Mathematical formulation of the problem
For the description of the ISD motion in the heliosphere we use a kinetic approach. In thisway we should calculate the ISD distribution function f d ( t, r , v ). The kinetic equation for f d ( t, r , v ) is: ∂f d ∂t + v · ∂f d ∂ r + F · ∂f d ∂ v = 0 , (1)where F is the sum of forces acting on the dust particles. On the right hand side of (1) wehave zero, because in the heliosphere one can neglect collisions between dust grains and 4 –their interaction with plasma protons and electrons (Gustafson 1994). In this article weconsider a stationary model of the solar magnetic field in the focusing phase, that is whythe solution of the kinetic equation is also stationary, ∂f d ∂t = 0: v · ∂f d ∂ r + F · ∂f d ∂ v = 0 . (2)The equation (2) requires boundary conditions to obtain a solution. In order to understandhow the ISD distribution is modified inside the region of the supersonic solar wind, weassume the ISD flow is undisturbed out of the Termination Shock (TS) - the shock wavethat restricts the region of the supersonic solar wind in the model of interaction betweensolar wind and interstellar medium. This helps to understand the modification of the ISDdistribution inside the TS as opposed to that in the heliospheric interface (Alexashov etal. 2016). We consider the TS as a sphere with radius r T S and formulate the boundarycondition as: f d ( r , v ) | r = r TS , v · e n > = f T S ( v ) , (3)where f T S ( v ) is the ISD distribution function on the TS and e n is the interior unit normalto the sphere. Below we discuss the form of function f T S ( v ) in more detail.To complete the correct mathematical formulation of the problem we should also set theboundary condition in the velocity space: f d ( r , v ) | v →∞ = 0 . (4)Note that each specific problem described by the formulation (2) - (4) is determined bya specific expression for the force term F ( r , v ) and for the boundary condition function f T S ( v ). Force analysis
Consider the Cartesian coordinate system as shown in Figure 1. Katushkina, Izmodenov(2019) provided the analysis of forces acting on the ISD particles. Four main forces acton the particles: the gravitational force F grav , the radiation pressure force F rad , the dragforce F drag due to interaction of dust grains with protons, electrons and neutrals and theelectromagnetic force F el . Estimates show (Gustafson 1994) that in the heliosphere wecan neglect the drag force.The expression for the gravitational force F grav is: F grav = − GM S r e r , (5)where G is the gravitational constant, M S is the mass of the Sun. 5 – Figure 1.
The coordinate system. The Sun is located at O, the velocity of the LISM v ISM is collinear to the Oy -axis. The Oz -axis coincides with the solar rotation axis. Sphericalcoordinates are introduced in the standard way.Since F grav is parallel F rad and both are proportional to r − , it is convenient to introducethe parameter β : β = | F rad || F grav | . (6)In this paper for the sake of simplicity we consider spherical particles. In this case β depends only on the star characteristics and particle mass m (see e.g. Katushkina, Iz-modenov 2019). Here we use the β = β ( m ) curve from Sterken et al. (2012) (green solidline in Figure 14). The resulting expression for the radiation pressure force is: F rad = β GM S r e r . (7)The magnetic field lines are frozen in the solar wind, that is why if we consider the referenceframe related to the solar wind, one can derive the expression for the electromagnetic forceusing relative (with respect to the solar wind) velocity v rel of dust particles: F el = qc m d ( v rel × B ) , (8)where v rel = v − v p is the dust particle velocity with respect to the solar wind, q is theparticle charge, c is the speed of light, m d is the dust grain mass, v p is the solar windvelocity, B is the solar magnetic field. The particle charge is expressed through the surfacepotential U d and the radius a : q = U d a , and we consider U d constant in the supersonicsolar wind (figure 2 from Alexashov et al. 2016, figure 2 from Slavin et al. 2012). Out of 6 –the region of the supersonic solar wind one should take account of the changes in valueof the surface potential, but it is beyond the scope of the present work. The dust grainmass m d = ρ d πa , where ρ d is the mass density of dust (here we consider astronomicalsilicates). We further assume uniform spherically symmetric solar wind: v p = v sw e r , andfor the solar magnetic field we use Parker’s model: B r = ± B E (cid:16) r E r (cid:17) , B φ = ∓ B E Ω r E v sw (cid:16) r E r (cid:17) sin θ, B θ = 0 , (9)where B E is the averaged solar magnetic field magnitude at the Earth’s orbit, r E is theastronomical unit, Ω is the angular velocity of the solar rotation. The sign ± denotes thechange in the polarity of the magnetic field across the heliospheric current sheet (HCS).Here for simplicity we assume a planar shape of the HCS (the plane Oxy in Figure 1).In reality there is a non-zero angle between the solar rotation axis and the magnetic axis,so the HCS has the ”ballerina skirt” shape. Here we also assume that the heliosphericmagnetic field is stationary: the HCS plane is at rest and in the region z > B r < , R ϕ > F ( r , v ) is: F ( r , v ) = ( β − GM s · e r r + qc m d (( v − v p ) × B ) . (10) Boundary condition
Let us assume that in the LISM there is a flux of ISD particles with the average velocity v ISM and dispersion of the v z velocity component. The interstellar magnetic field hasspatial and temporal inhomogeneities which act as sources for the acceleration of ISDparticles (Hoang et al. 2012). This acceleration is the reason for variations in the velocitiesof individual dust particles and, therefore, appearance of dispersion in the ISD velocitydistribution. Below we demonstrate that relatively small values of dispersion of the v z component significantly influence the results. Then the expression for f T S ( v ) is: f T S ( v ) = n ISM δ ( v x ) δ ( v y + v ISM ) 1 σ z √ π exp (cid:18) − v z σ z (cid:19) , (11)where δ is the Dirac delta-function, σ z is the dispersion of the v z component. As σ z → v ISM : f T S ( v ) = n ISM δ ( v x ) δ ( v y + v ISM ) δ ( v z ) , (12) 7 –and the formulation of the problem is identical to the one given by Mishchenko et al.(2020).Slavin et al. (2012) take account of the dispersion by addition of the supplementaryvelocity component lying in the plane perpendicular to the direction of the interstellarmagnetic field. This supplementary velocity component has constant absolute value (3km/s) and random direction in the above-mentioned plane. In the present work we modelthe dispersion of the v z component using the normal distribution. Dimensionless formulation of the problem
As a characteristic distance we consider L = GM S v ISM and as a characteristic velocity – v ISM .Since the problem is linear and homogeneous in f d ( r , v ) we can substitute f d → f d n ISM andeliminate n ISM in (11). The dimensionless formulation of the problem (2) - (4), (10), (11)is: ˆ v · ∂ ˆ f d ∂ ˆ r + ˆ F · ∂ ˆ f d ∂ ˆ v = 0 , ˆ f d (ˆ r , ˆ v ) | ˆ r =ˆ r TS , ˆ v · e n > = δ (ˆ v x ) δ (ˆ v y + 1) σ z √ π exp (cid:16) − ˆ v z σ z (cid:17) , ˆ f d (ˆ r , ˆ v ) | ˆ v →∞ = 0 , (13)where ˆ r = r L , ˆ v = v v ISM , ˆ f d = f d v ISM , ˆ F = F L v ISM , ˆ σ z = σ z v ISM , ˆ r T S = r TS L . The expression forthe sum of forces (10) in the dimensionless form is:ˆ F = ( β − e r ˆ r + sgn (ˆ z ) v em v ISM (cid:18) v ISM v sw ˆ v − e r (cid:19) × (cid:18) − L Ω L e r ˆ r + sinθ ˆ r e ϕ (cid:19) , (14)with L Ω = v sw Ω , v em = qB E Ω r E c m d . The dimensionless formulation of the problem contains fivedimensionless parameters:ˆ σ z , β, ε = v em v ISM = 3 U d B E R E Ω4 πc ρ d a v ISM , v
ISM v sw , L Ω L . (15) Trajectories in the plane of symmetry
Let us consider the projection of the force ˆ F on the x -axis:ˆ F x = ( β −
1) ˆ x ˆ r + sgn (ˆ z ) ε ˆ x ˆ z ˆ r + sgn (ˆ z ) ε v ISM v sw (cid:18) − L Ω L ˆ v y ˆ z − ˆ v z ˆ y ˆ r − ˆ v z ˆ x ˆ r (cid:19) . (16)Since v ISM v sw ≈ . <<
1, one can neglect the corresponding term, and the expression forˆ F x in a simplified form is: 8 – Figure 2.
The computational domain is a square with the side 2ˆ r T S = r TS L which isdivided into rectangular cells ∆ˆ y × ∆ˆ z .ˆ F x = ( β −
1) ˆ x ˆ r + sqn (ˆ z ) ε ˆ x ˆ z ˆ r . (17)Therefore, dust particles with initial parameters:ˆ x | T S = 0 , ˆ v x | T S = 0 , (18)can’t leave the plane ˆ x = 0 under the action of the force (17) according to Picard’s existenceand uniqueness theorem. For simplicity we consider only such trajectories in this article. Monte-Carlo approach
To solve the kinetic equation we use the Monte-Carlo method. The computational domainis divided into rectangular cells ∆ˆ y × ∆ˆ z (Figure 2) and ∆ˆ z << ∆ˆ y because the singularlayers found in Mishchenko et al. (2020) are oriented horizontally. Moreover, since theirthickness approaches zero one should decrease ∆ˆ z in order to detect these peculiarities bythe Monte-Carlo modeling.For a dust particle we generate randomly its initial velocity and location on the sphere withradius r = r T S according to the distribution function f T S ( v ) from (13). During the motionof the particle in the heliosphere we record the time t i of the particle in the computationaldomain cells ( t i = 0 if particle does not cross the corresponding cell). Then, by thedefinition of the distribution function and number density, and the law of large numberswe have: 9 –ˆ F N N (cid:88) i =1 t i ∆ˆ r c ∆ˆ v c → f d (ˆ r c , ˆ v c ) (19)ˆ F N N (cid:88) i =1 t i ∆ˆ r c → n d (ˆ r c ) , (20)where N is the number of particles, ∆ r c ∆ v c is the cell volume in the phase space, ˆ F is theflux of the dust particles through the outer surface per unit of time in the dimensionlessform: ˆ F = π (cid:90) − π (cid:90) (ˆ v · e n ) > (ˆ v · e n ) ˆ f T S (ˆ v ) d ˆ v ˆ r T S dϕ = 2ˆ r T S (21)
Technical characteristics
In this paper we consider particles with the radius a = 0 . µm . For these particles β = 1and consequently the gravitational and radiation pressure forces cancel out in (14).For computations we use the following values of the parameters: r T S = 100 a.u, v ISM =26 . M S = 2 · kg, v sw = 400 km/s, Ω = 2 . · − U d = +3 V, B E = 30 µ G, R E = 1 a.u., ρ d = 2500 kg/m .For all figures with results in this paper, unless otherwise specified, the cell size inside thecomputational domain is 0 . a.u. × . a.u. in the Oy - and Oz -directions, respectively.To solve the system of ODEs for the trajectory of a particle the fourth order Runge-Kuttamethod was used.Note that the selected fixed location of the HCS corresponds to the case when all ISDparticles are attracted to the HCS (focusing phase). In order to understand it let usconsider the z -axis projection of (14):ˆ F z = ( β −
1) ˆ z ˆ r + sgn (ˆ z ) ε (cid:18) − ˆ x + ˆ y ˆ r + v ISM v sw (cid:18) − L Ω L ˆ v x ˆ y − ˆ v y ˆ x ˆ r + ˆ v x ˆ x + ˆ v y ˆ y ˆ r (cid:19)(cid:19) , (22)where again v ISM v sw ≈ . <<
1, that is why at large heliospheric distances the leadingterm is: − sgn (ˆ z ) ε ˆ x + ˆ y ˆ r . (23)It is seen, that in the case ˆ z > z -axis component ˆ F z <
0, so the ISD particles areattracted to the HCS. In the case ˆ z < F z > Figure 3.
Map of the density distribution in the plane X = 0 in the case withoutdispersion. Yellow color envelopes correspond to caustics. Relative statistical error islimited by 2-3 % at each point. The number of trajectories N = 2000000. For the sake ofcomparison, the panel at the right bottom presents the results of Mishchenko et al. (2020)obtained for the same conditions. The radius of particles is 0.37 µ m. RESULTS
Singularities in density
It was shown by Mishchenko et al. (2020) that in the case of zero dispersion the ISDtrajectories form caustics at which the number density of ISD is infinite. A caustic is anenvelope of the ISD trajectories. By definition, every segment of a caustic is tangent to aninfinite number of the ISD trajectories, that is the reason for density singularities origin.The distribution of the dust density has multiple singularities. This result was obtained bythe Lagrangian approach. In this Section we demonstrate that the singularities can be alsoobtained by the Monte-Carlo approach (although, perhaps, at the cost of computationalefficiency).Figure 3 shows the map of the ISD density distribution as well as the ISD streamlines.The map shows the region in the vicinity of the HCS. Symmetrical yellow lines in thisFigure are the above-mentioned caustics. In the case of the Monte-Carlo simulation theyrepresent the thin regions where a sharp density peak is found (Figure 4). Inside the areadelimited by the caustics there is a complex structure of the ISD density distribution withmany local peaks.Since the computational domain consists of finite size cells, high spatial resolution of thenumerical grid is required to detect the density singularities with high precision by theMonte-Carlo modeling. Figure 5 shows how the ISD density distribution along the line 11 –
Figure 4.
The distribution of ISD density along the line ( X = 0 , Y = 2). The cell sizeis 0 . × .
001 a.u. Relative statistical error is limited by 2-3 % at each point. Thenumber of trajectories N = 2000000. The radius of particles is 0.37 µ m. Figure 5.
The distribution of ISD density along the same line as in Figure 4. Thefour lines of the different colour correspond to the different cell sizes in the z -direction(∆ y = 0 . z = { .
05 a.u., 0 .
01 a.u., 0 .
005 a.u., 0 .
001 a.u. } ). Relative statisticalerror is limited by 2-3 % at each point. The number of trajectories N = 2000000. Theradius of particles is 0.37 µ m. 12 – Figure 6.
The tube of trajectories from a small region on the outer boundary. The tubeis compressed up to a thousand times at small heliocentric distances. The point of minimalwidth corresponds to a point on the caustic. Small boxes demonstrate trajectories view ata large scale. The box sizes are 0 . × .
025 a.u. The radius of particles is 0.37 µ m.( X = 0 , Y = 2) changes with variation of cell size ∆ z . The ISD density at cells containingcaustic points increases with decreasing ∆ z and, therefore, the ISD density singularitiesare located at these cells.A simple explanation for the formation of the caustics is as follows. Let us look at theparticles originating from a small region on the boundary of the computational domain(Figure 6). This flux tube is compressed with decreasing heliocentric distance and reachesits minimal width (approaching zero) exactly at the caustic points. Considering the con-servation of mass for a flux tube: n v Σ = n v Σ a minimal value of tube width Σ corresponds to the maximal value of the density n ,because the value of the y -component velocity v is approximately constant. Effects of velocity dispersion
In this paper we mainly consider the dispersion of the v z component. We do not considerthe dispersion of the v x component at all because we only study the plane of symmetry X = 0, and Figure 11 shows that the dispersion of the v y component has less impact onthe density distribution than the dispersion of the v z component.To explore the effect of velocity dispersion we performed the calculations for the set of ˆ σ z values: 0 , . , . , .
04. Figure 7 presents the density maps obtained for the four values ofˆ σ z . With increasing ˆ σ z the density maxima are smeared and their singularities disappear. 13 – Figure 7.
Comparison of the maps of the ISD density distributions with different disper-sion ˆ σ z : 0 , . , . , .
04. With increasing dispersion ˆ σ z the caustics are smeared and thedensity singularities disappear. Relative statistical error is limited by 2-3 % at each point.The number of trajectories N = 2000000. The radius of particles is 0.37 µ m. Figure 8.
Comparison of the ISD density distributions along the line ( X = 0 , Y = 2)with different dispersion ˆ σ z : 0 , . , . , .
04. The structure of the density distributionchanges drastically with variation of dispersion. Relative statistical error is limited by 2-3% at each point. Th number of trajectories N = 2000000. The radius of particles is 0.37 µ m. 14 – Figure 9.
Trajectories of particles originating from a small region on the boundary ofthe computational domain for different values of ˆ σ z . Particle trajectories are scattered fornon-zero dispersion, and no singularities appear. The radius of particles is 0.37 µ m. Figure 10.
Comparison of the fluxes carried by the ISD trajectories originating from asmall region on the boundary of the computational domain with different dispersion ˆ σ z :0 , . , . , . , .
04. For each computation we have the convergence in cell size,Runge-Kutta integration steps and the number of simulated particles. The maximal valueof fluxes is approximately inversely proportional to the dispersion ˆ σ z . For computationwithout dispersion the cell size in the z -axis direction is 10 − (a.u.). Relative statisticalerror is limited by 5 % at each point. The radius of particles is 0.37 µ m. 15 – Figure 11.
Comparison of density map distributions with different dispersion of the v y velocity component: 0 , . , . , .
04. The shape of the overdensity regions remainsvirtually unchanged with the variation in the dispersion. Relative statistical error is limitedby 5-6 % at each point. Number of trajectories N = 200000. Radius of particles is 0.37 µ m.The regions of overdensity remain only in the vicinity of the HCS. This clumping up ofdust particles is associated with an increase in the magnitude of the Lorentz force at smallheliocentric distances, which leads to decrease in the amplitude of the particle oscillationsaround the HCS. Then, the gross tendency is for the ISD to converge to the HCS planeand so the regions of overdensity appear. In Figure 8 we can see how the density atcells containing caustic points changes quantitatively with variation of ˆ σ z . Small values ofdispersion in the boundary velocity distribution drastically change the density distributionin the heliosphere.The reason for the disappearance of the singularities is clearly seen from Figures 9, 10. Fig-ure 9 shows trajectories of particles originating from a small region on the outer boundaryof the computational domain for different ˆ σ z values. As it was mentioned above, singu-larities appear where the width of the flux tube approaches zero. Particle trajectories arescattered for non-zero dispersion, and therefore no singularities appear. Density flux dis-tributions in the vicinity of the point corresponding to the caustic by particles originatingfrom a small region on the outer boundary are demonstrated in Figure 10 for different ˆ σ z values. One can appreciate that the maxima of these density flux distributions are ap-proximately inversely proportional to the value of the dispersion ˆ σ z . Thus even extremelysmall values of dispersion drastically influence the ISD distribution inside the heliosphere.In order to study the influence of the dispersion of the v y component on the densitydistribution instead of the expression (11) we should use the following boundary conditionfunction: 16 – f T S ( v ) = n ISM δ ( v x ) 1 σ y √ π exp (cid:18) − ( v y + v ISM ) σ y (cid:19) δ ( v z ) , (24)which in the dimensionless form is: f T S (ˆ v ) = δ (ˆ v x ) 1ˆ σ y √ π exp (cid:18) − (ˆ v y + 1) σ y (cid:19) δ (ˆ v z ) . (25)Figure 11 shows the comparison of density distributions for the cases with different values ofdispersion ˆ σ y of the v y component (ˆ σ y : 0 , . , . , . v z component, one can concludethat the dispersion of the v z component has greater impact on the density distribution thanthe dispersion of the v y velocity component. This is because the regions of overdensity arestretched along the Oy -axis and, therefore, in the stationary case small variations in the v y component can’t significantly influence the ISD density distribution. CONCLUSION
In this paper we demonstrated that the singularities of the ISD density in the heliosphere,discovered using the Lagrangian approach in Mishchenko et al. (2020), can also be foundby the Monte-Carlo simulations. This requires super-small computational cells. In ourcalculations the required size of a cell (in the z -direction ) is 10 − a.u. Having a suchsize of the cells in the whole domain is computationally unrealistic. Weaker resolution (i.e.larger cells) does not allow to find the caustics.Dispersion was introduced as a normal distribution of one of the velocity component.It was shown that the density singularities are smeared due to dispersion. The regionsof overdensity are smoothed and remain only in the vicinity of the heliospheric currentsheet. It is known (Hoang et al. 2012) that the velocity dispersion can reach valuesof approximately 15 % due to spatial and temporal inhomogeneities in the interstellarmagnetic field. Significant qualitative and quantitative changes in the density distributionemerge even for 5 % dispersion as it was shown. Thus, the velocity dispersion is anextremely important effect that strongly influences the ISD density distribution inside theheliosphere.In the future we plan to develop our model to the case of the time-dependent solar magneticfield in accordance with the 22-year solar cycle (in this paper we considered the solarmagnetic field just in one focusing phase (Mann 2010)). Certainly this is a highly importanteffect that has a major impact on the ISD density inside the heliosphere and which isnecessary to take into consideration. Acknowledgements
17 –The authors are grateful to the Government of Russian Federation and the Ministry of Sci-ence and Higher Education for the support by grant 075-15-2020-780 (N13.1902.21.0039).We thank D. B. Alexashov, I. Baliukin and A. Granovskiy for useful discussions and forthe help with preparation of the manuscript. This work is supported by grant 18-1-1-22-1of the ”Basis” Foundation. 18 –
REFERENCES Alexashov et al. (D. B. Alexashov, O. A. Katushkina , V. V. Izmodenov, P. S. Akaev),MNRAS , 2553 (2016).2.
Altobelli et al. ( N. Altobelli, S. Kempf, H. Kr¨uger, M. Landgraf, M. Roy, E. Gr¨un),Journal of Geophysical Research , 7102 (2005).3.
Altobelli et al. (N. Altobelli, V. Dikarev, S. Kempf, R. Srama, S. Helfert, G. Moragas-Klostermeyer, M. Roy, E. Gr¨un), Journal of Geophysical Research , 7105 (2007).4.
Baliukin et al. ( I. I. Baliukin, V. V. Izmodenov, E. M¨obius, D. B. Alexashov, O. A.Katushkina, H. Kucharek), Astrophys. J. , 119 (2017).5.
Bertaux, Blamont ( J. L. Bertaux, J. E. Blamont), Astron. Astrophys. , 200 (1971).6. Bertaux, Blamont ( J. L. Bertaux, J. E. Blamont), Nature , 263 (1976).7.
Witte et al. (M. Witte, H. Rosenbauer, E. Keppler, H. Fahr, P. Hemmerich, H.Lauche, A. Loidl, R. Zwick), Astron. Astrophys. , 333 (1992).8. Witte (M. Witte), Astron. Astrophys. , 835 (2004).9.
Gr¨un et al. (E. Gr¨un, B. Gustafson, I. Mann, M. Baguhl, G. E. Morfull, P. Staubach,A. Taylor, H. A. Zook), Astron. Astrophys. , 915 (1994).10.
Gustafson (B. A. S. Gustafson), Ann. Rev , 553 (1994).11. Draine (B. T. Draine), Space Sci. Rev. , 333 (2009).12.
Zirnstein et al. (E. J. Zirnstein, J. Heerikhuisen, H. O. Funsten, G. Livadiotis, D. J.McComas, N. V. Pogorelov) Astrophysical Journal Letters , 30 (2016).13.
Izmodenov et al. (V. V. Izmodenov, Y. G. Malama, A. P. Kalinin, M. Gruntman, R.Lallement, I. P. Rodionova), Astrophys. Space Sci. , 71 (2000).14.
Izmodenov, Alexashov (V. V. Izmodenov, D. B. Alexashov) Astrophys. J. Suppl. Ser., , (2015).15.
Izmodenov, Alexashov (V. V. Izmodenov, D. B. Alexashov) Astron. Astrophys., , (2020).16.
Katushkina et al. (O. A. Katushkina, V. V. Izmodenov, D. B. Alexashov, N. A.Schwadron, D. J. McComas), Astrophys. J. Suppl. Ser., , (2015).17.
Katushkina, Izmodenov (O. A. Katushkina, V. V. Izmodenov), MNRAS , 4947(2019).18.
Qu´emerais et al. (E. Qu´emerais, B. R. Sandel, V. V. Izmodenov, G. R. Gladstone),
Cross-Calibration of Far UV Spectra of Solar System Objects and the Heliosphere (2013).19.
Landgraf et al. (M. Landgraf, W. J. Baggaley, E. Gr¨un, H. Kr¨uger, G. Linkert),Journal of Geophysical Research , 10343 (2000).20.
Landgraf et al. (M. Landgraf, H. Kr¨uger, N. Altobelli, E. Gr¨un), Journal of Geophys-ical Research , 8030 (2003). 19 –21.
Levy, Jokipii (E. H. Levy, J. R. Jokipii), Nature , 423 (1976).22.
McComas et al. (D. J. McComas, M. Bzowski, P. Frisch, S. A. Fuselier, M. A. Kubiak,H. Kucharek, T. Leonard, E. M¨obius et al.), Astrophys. J. , 28 (2015).23.
Mann (I. Mann), Annu. Rev. Astron. Astrophys , 173 (2010).24. Mathis et al. (J. S. Mathis, W. Rumpl, K. H. Nordsieck), Astrophys. J. , 425(1977).25.
Moebius et al. (E. Moebius, P. Bochsler, M. Bzowski, G. B. Crew, H. O. Funsten, S.A. Fuselier, A. Ghielmetti, D. Heirtzler et al.), Science , 969 (2009).26.
Mishchenko et al. (A. V. Mishchenko, E. A. Godenko, V. V. Izmodenov), MNRAS , 2808 (2020).27.
Osiptsov (A. N. Osiptsov), Astrophysics and Space Science , 377 (2000).28.
Pogorelov et al. (N. V. Pogorelov, J. Heerikhuisen, G. P. Zank, S. N. Borovikov, P.C. Frisch, D. J. McComas), Astrophys. J. , 104 (2011).29.
Slavin et al. (J.D. Slavin, P.C. Frisch, H.-R. M¨uller, J. Heerikhuisen, N. V. Pogorelov,W. T. Reach, G. P. Zank), Astrophys. J. , 46 (2012).30.
Sterken et al. (V. J. Sterken, N. Altobelli, S. Kempf, G. Schwehm, R. Srama, E.Gr¨un), Astron. Astrophys. , A102 (2012).31.
Sterken et al. (V. J. Sterken, A. J. Westphal, N. Altobelli, D. Malaspina, F. Postberg),Space Sci. Rev. , 7, 43 (2019).32.
Strub et al. (P. Strub, H. Kr¨uger, V. J. Sterken), Astrophys. J. , 140 (2015).33.
Strub et al. (P. Strub, V. J. Sterken, R. Soja, H. Kr¨uger, E. Gr¨un, R. Srama), Astron.Astrophys. , A54 (2019).34.
Hoang et al. (T. Hoang, A. Lazarian, R. Schlickeiser), Astrophys. J. , 54 (2012).35.
Czechowski, Mann (A. Czechowski, I. Mann), Astron. Astrophys.410