Calculation of a magnetic force acting on small superconducting celestial bodies

Recent discoveries of superconducting phases in the samples of meteorites suggest the possibility of a natural occurrence of superconducting state in space. Superconductors are known to exhibit interesting behaviours when subjected to external magnetic fields, such as levitation. Similar force may act on a superconducting bit in space. The goal of this paper is to quantify this force and assess its effects. Several scenarios in which a superconducting bit can be produced and interact with a magnetic field in space are suggested. The force acting on a superconductor in different conditions is calculated with numerical simulations. The dependence on a magnetic flux density, its gradient, and the geometry and the properties of the superconductor are found. The empirical formulas are derived and used to calculate a magnetic force. The resultant force is extremely weak in all analysed scenarios. It is found that its strength decreases rapidly with the distance from the source of the magnetic flux. Its effect on trajectory of the superconductor is almost negligible. Some possibilities of increasing its strength and the effects are considered.


Introduction
Naturally occurring superconducting phases were recently found by Wampler et al. [1] in the samples from two meteorites-iron meteorite Mundrabilla and ureilite GRA 95205 using magnetic field modulated microwave spectroscopy [2]. Some other hints were observed earlier by Guenon et al. [3]. This discovery prompts to analyse the possible effect of superconductivity on the behaviour of celestial bodies and quantify the strength of a magnetic force acting on possible superconducting bodies in space.
The presence of natural superconductors on Earth is rare, and the search for them yielded mixed results [4]. Superconductivity in natural terrestrial minerals was observed in covellite by di Benedetto et al. [5] cooled down below its T c of 1.63 K, and in miassite, with T c of 5.4 K [6]. Some elements naturally present in native form are superconductors too. Notably, T c of tin is 3.72 K and of lead is 7.19 K [7]. It was the alloy of these elements (along with indium), which was detected by Wampler et al. [1]. Another candidate for a space superconductor are some types of diamonds. They are quite abundant, both in stardust with sizes of approximately 2 nm [8] and as larger grains, up to 100 μm [9]. Boron-doped diamonds a e-mail: lukasz.tomkow@cern.ch (corresponding author) an example of a planet ejected from a planetary system and maintaining bodies of water on its surface. Due to the thick atmosphere, such a hypothetical planet would have an effective temperature of approximately 30 K, giving equilibrium temperature around it well below T c of MgB 2 . Meanwhile, high temperature on the surface would be kept thanks to dense and thick atmosphere, nuclear reactions in its interior and residual heat from the time of its formation. The latter implies the presence of the convective currents, which could lead to the generation of a magnetic flux. Another possibility of heating the interior are tidal forces possibly acting on moons of unbound planets [39]. Brown dwarfs also can act as unbound sources of strong magnetic fluxes (with densities up 4.1×10 −1 T) with relatively low surface temperatures estimated as slightly above 800 K [40].
The goal of this paper is to analyse the behaviour of a small superconducting bit in a magnetic flux that could be generated by a larger celestial body. In order to achieve this goal, a number of numerical simulations are performed to develop simple empirical formulas for a magnetic force acting on a superconductor interacting with a magnetic flux with a given density and its gradient. The formula should also take into consideration the properties of a superconducting material, such as critical current density J c and T c , its geometry. For continuity, the results from numerical calculations used to find the formulas for magnetic force are presented along its formulation in methods. Then, the magnetic force acting on a superconducting bit is calculated in different scenarios. The analysed scenarios include orbiting around a Sun-like star, a magnetar, a brown dwarf, and a rogue planet exhibiting an Earth-like magnetosphere.

Methods
To calculate a magnetic force acting on a superconductor in different conditions, a numerical model is developed in COMSOL Multiphysics. Magnetic field formulation package is employed, which uses H-formulation for calculation of the magnetic fields and current densities, as described further. Several assumptions are made to simplify the numerical model and ease the extraction of dependencies on several factors. The calculations are performed in 2D with circular and rectangular geometries presented schematically in Fig. 1. In the actual model, the environment region is much larger than a superconducting region. It is modelled as a vacuum with relative magnetic permeability μ r equal to 1. Because of that, the magnetic induction B on the edges of the modelled region can be found as μ 0 H , where μ 0 is magnetic permeability of vacuum and H is magnetic field.
Prescribed value of B is applied as boundary condition on the edges. It is assumed that it has only one nonzero component B y . Analyses are performed for different combinations of background magnetic induction B b and its gradient along axis x, ∇ B y(x) The example value of B y versus time is presented in Fig. 2a. Initially B y is 0 on both edges. During the first 0.25 s, it is increased to the value B b used at a given modelling run, resulting in uniform applied B y across all boundaries. It is maintained as such for 0.75 s to allow for relaxation. Then, ∇ B y(x) is gradually increased. After each increase, the applied B y is constant for the period of 1 s and the total value of magnetic force is read at the end of this time. ∇ B y(x) is changed to have positive and negative values as shown in Fig. 2b to account for the possible hysteretic effects.  H-formulation is based on the work by Pecher et al. [41], which in turn is based on Maxwell's formulas. H-formulation is widely used in modelling of superconductors [42], and it has been successfully applied to model the levitation behaviour of a superconductor by Sass et al. [43]. While applying the formulation, it is assumed that the excitation frequencies are low, allowing to use the quasi-static electromagnetic model. With these assumptions, Faraday's equation takes form of formula 1. The magnetic field is denoted as H, the electric field is E, t is time, μ 0 is magnetic permeability of vacuum and J is electric current density.
To satisfy Ampere's law, J is expressed with formula 2, also being a quasi-static approximation.
Two different formulas are used to find E in different regions of a model, as shown in Fig. 1. In a non-superconducting environment region, representing vacuum surrounding a superconductor, E is found with formula 3. σ here is electrical conductivity. In H-formulation, slightly unphysically, a very small value (less than 1000 S/m) of σ has to be attributed to the environment region to maintain a numerical stability.
Current-voltage characteristics of a superconductor is highly nonlinear. To account for that in a superconducting region, E is calculated depending on a current density J with the power law (formula 4). In this formula, J c is a critical current density. During the calculations of the magnetic force, it is assumed as constant for a given case, to ease the extraction of the effect of its value on the force. In further calculations, the dependence on the magnetic induction and temperature is introduced. Similarly, the exponent of power law n is assumed to have a constant value of 40. Threshold electric field E 0 is assumed as 100 μV m −1 [44]. To account for the finite size of the superconducting bit, a constraint setting the total value of the electric current in the superconducting region to 0 is imposed. The numerical mesh used in the model is very dense in the region close to the surface of superconductor, where the shielding electric currents are expected to appear.
Lorentz force density f L is calculated using formula 5. Since a 2D model is used, only a total force per unit length F m /L can be found. This value is the integral of f L over the surface of superconducting region. In order to find the total value of the magnetic force F m , F m /L has to be multiplied by L, taking into consideration the shape changes along z axis, perpendicular to the x y plane shown in Fig. 1.
The origin of net F m acting on a superconductor is the inhomogeneity of the density of electric current, which arises from the inhomogeneity of B b . Thus, the net F m appears only in a non-uniform background magnetic flux. The strength of F m is affected by several factors. The ones considered in this paper are strength and degree of non-uniformity of B b , material properties of the superconductor represented by J c , and the geometry of the superconducting bit. To quantify their effect, the numerical computations are performed in the wide range of possible values and the resulting F m is calculated. The base values in the computations are: background magnetic flux density of 1 mT, diameter of superconductor (for circular geometry) or side length (for rectangular geometry) of 0.02 m and critical current density J c of 2 × 10 8 A/m 2 .
Dependence of F m on B b and its gradient is presented in Fig. 3. For a given B b , the strength of F m rises linearly with its gradient, as shown in Fig. 3a. No significant dependence of strength and direction of F m on an angle between B b and ∇B b is observed. It is important in the considerations of the effect of geometry-different formulas apply for the dimensions in different orientations, as shown further. The net F m always acts towards decreasing B b . In further analysis, the values of F m obtained for a given parameter are collectively considered in terms of a linear coefficient a to help in the representation of dependencies. a of F m L (∇ B) increases linearly with B b in a wide range of its strength.
Some deviations from the linear dependence of a on B b can be seen in Fig. 3b for very low and high values of B b . The source of the deviation in the lower range of B b is the reversal of direction of B y occurring in the superconductor region even a when a very small gradient is applied. Additionally, in such conditions the shielding currents flow in the layer much thinner than the resolution of a numerical mesh. Thus, it can be safely assumed that the validity of linear fit extends into the low B b range. In the case of high B b , a magnetic flux can penetrate the superconductor beyond the surface layer, leading to more complex patterns of the shielding currents. The values of B b , in which the deviation is significant, are higher than the ones realistically encountered in space by a superconductor. Thus, the linear fit is sufficient for the purpose of the analysis carried out in this paper. Hysteretic effects are negligible in every analysed case.
Values of |F m |/L versus critical current density J c are presented in Fig. 4. The effect of changing J c even by several orders of magnitude appears to be rather weak in Fig. 4a, but the analysis of the directional coefficients reveals a clear dependence shown in Fig. 4b. A good fit is obtained with formula 6. The values of fitting parameters for b 1 to b 4 , respectively, are: 0.248, 7.95 × 10 −9 , -0.075 and 1.48 × 10 −2 for J c given in A/m 2 . The observed dependence comes from lower local J and the deeper penetration of a superconductor by the magnetic flux, resulting in the lower local values of f L acting on superconductor and, ultimately, a smaller net F m . It can be expected that with J c approaching 0, F m will disappear whatsoever, as there will be no shielding currents.
The effect of the superconductor geometry on F m is analysed for two geometries-circular and rectangular. Dependence of F m on radius r of a circular superconductor is presented in  Figure 5a shows the results of numerical calculations for different values of r . The value of directional coefficient increases with the square of r , as can be seen in Fig. 5b. In a magnetic field with a constant gradient, the total difference of magnetic flux density across the superconductor rises linearly with increasing r . Similarly, the circumference of the ring in which the shielding currents are flowing also increases linearly, as it is proportional to radius of a circle. The combined effect of these factors gives rise to the observed quadratic dependence.
In the case of the rectangular geometry, the increase of the dimension of the superconductor parallel to ∇B b has a very similar effect, as shown in Fig. 6. In the presented analysis, this dimension is the width of the sample and is represented as ratio k w = w/ h. As in previous case, the increasing difference of B b across the superconductor, and of the surface where the non-symmetric currents appear, combine to form a quadratic dependence (Fig. 6b). The effect of changing the height of the sample (the dimension perpendicular to ∇B b ) can be described by a logarithmic function defined with formula 7. The results of fitting and calculations for changing height are shown in Fig. 7. The fitted parameters c 1 , c 2 and c 3 are 0.106, 21.57 and 11.78, respectively.
Based on the above considerations, formula 8 is used to find |F m |/L acting on a circular crosssection of a superconductor. Analogically, formula 9 is used for rectangular cross-sections. The value of constant C c is approximately −2.49 × 10 6 and C r is −2.67 × 10 5 , when all parameters are given in SI units. As mentioned, to obtain the formulas for force acting on 3D objects, formulas 8 and 9 have to be integrated along z axis. In the case of a rectangular prism, the logarithmic dependence has to be considered as well. The resultant formulas are 10 for a sphere and 11 for a rectangular prism. The presented formulas are valid for |B b | lower than approximately 1 T. They do not consider the effects of field cooling of a superconductor, and the fitted values of parameters are only approximations. Nevertheless, the formulas are sufficiently accurate to provide the realistic values of forces acting on a superconductor in space and to allow the prediction of its behaviour.
In the view of the reasoning presented in Introduction, the most realistic candidate for a superconductor naturally present in space is MgB 2 and its material properties will be used in further calculations. Since it is characterised by low anisotropy [45], J c is assumed to depend only on magnetic flux density B and temperature T . The values are interpolated based on experimental results gathered by Buzea and Yamashita [46] and shown in Fig. 8. The analysed superconductor is assumed to be a spherical piece of pure MgB 2 with the radius of 0.02 m. A temperature of a superconducting bit is assumed to be equal to equilibrium temperature T eq at a given location, found with formula 12. In this formula, T body is an average surface temperature of a celestial body interacting with a superconductor, which is treated as a black body. R is the celestial body's radius, s is the distance between the interacting objects and A B is Bond albedo of a superconductor, assumed as 0.25. The values of T body are listed in Table 1.
Calculation of B is based on a simple numerical model, assuming that a celestial body is a spherical region with an electric current flowing around its axis in angular direction. In other words, the body is treated as a spherical electromagnet. Current density is assumed to be uniform across the volume of the body and its value is selected so the maximum |B| on the surface matches the assumed one, according to Table 1. Further calculations are performed only on the surface perpendicular to the axis of the current, where only a single component B z of B exists. The obtained values of magnetic flux density are then fitted to find the exact value at a given distance from a body. Gradient is calculated numerically based on the interpolated data. Table 1 lists several feasible scenarios in which a superconducting bit can interact with bodies generating magnetic flux. Scenario A assumes that a superconductor is interacting with a Sun-like star with the magnetic flux on its surface chosen as the average of a typical range. In scenario B, an interaction with a magnetar is considered, including extremely high magnetic flux densities and temperatures. In scenario C, a relatively low-temperature brown dwarf is selected as a magnetic flux source. Finally, in scenario D a rogue planet with surface magnetic flux density similar to Earth is considered, the major difference being significantly lower effective surface temperature than Earth's.
In each scenario, the values of F m obtained with formula 10 are compared with gravitational force F G . Calculation of F G is performed using a standard formula (13). G is gravitational constant, d is distance between centre of mass of an orbited body and a superconductor. The ratio |F m | / |F G | is calculated at each respective point. At all considered scenarios, the two  Fig. 9 a Equilibrium temperature around a rogue planet, b resultant critical current density forces are acting in opposite directions-F G pulls a superconducting bit towards a larger body, while F m repels it.

Results and discussion
The calculated values of F m acting on a superconductor are very small at each scenario. The strongest F m is found with the assumptions from scenario D (rogue planet). The values obtained with this scenario are presented in the following figures. Equilibrium temperature and critical current density are shown in Fig. 11. In all scenarios, the strength of B has a negligible effect on critical current density, and its value is almost completely determined by temperature. Critical temperature of bulk MgB 2 at the conditions of almost zero magnetic field is 39 K. However, according to the data from Buzea and Yamashita [46], the temperature at which J c of MgB 2 reaches the values higher than 5 × 10 6 A/m 2 is 31.05 K. The considerations shown in Fig. 4 indicate that when J c is below this value, there is almost no magnetic force present. Assuming a uniform temperature distribution of a rogue planet surface, a temperature of 31.05 K is achieved at the distance of approximately 7340 km from its centre, approximately 340 km from its surface (Fig. 9a). Critical current density increases with the distance and decreasing equilibrium temperature, as shown in Fig. 9b.
B and ∇B quickly decrease with the distance from the planet centre-magnetic induction decreases as 1/d 3 Fig. 11 a Strength of force due to magnetic interaction between superconductor and a magnetic field of a rogue planet, b ratio between the strength of magnetic force and gravity in the considered region is 5 × 10 −5 T (Fig. 10a) and maximum ∇B is 1.05 × 10 −11 T/m (Fig. 10b). It leads to extremely small values of F m .
As shown in Fig. 11a, despite the relatively low critical current density, F m is strongest in the region closest to the planet, almost immediately after the equilibrium temperature is below the critical temperature. The maximum found value of F m is only 4 × 10 −15 N for the considered sphere of MgB 2 with r of 0.02 m. With the used assumptions, both the gravitational force and the magnetic interaction force depend on r as r 3 . Thus, the ratio between the two will remain the same regardless of the size of a superconducting bit. Assuming that the superconductor is initially on a circular orbit in the region where F m is maximum, the operation of the force will increase the circular orbit radius by 50 nm, making it almost negligible. In the considered geometry, an acceleration due to F m is dependent only on the position of a superconductor and not its mass, giving the maximum value of 4.62×10 −14 m/s 2 The results of calculations for all scenarios are gathered in Table 2. The obtained values of the magnetic interaction force are even lower for other bodies. With the force decreasing as 1/d 7 , the larger distances of superconducting transition present in other scenarios mean that even the extremely strong sources of magnetic flux such as magnetars fail to exert a meaningful force on a superconductor. On the other hand, such dependence opens some possibilities to increase the effect of superconductivity on a superconductor trajectory. One of the assumptions taken during the calculations of F m is that the temperature of a superconductor is equal to an equilibrium temperature at a given location. It is not strictly true, as each body possess a certain thermal inertia. Therefore, it is possible for a sufficiently large and fast-moving superconductor passing by an object to maintain a superconducting state in the region with higher magnetic flux density and thus to experience a stronger force. In the considered example of a rogue planet, the closest possible approach to d equal to an assumed atmosphere radius leads to a twofold increase of F m when compared with the case not considering heat capacity of a superconductor.
The assumed uniform, spherical form of a superconductor is rather unlikely. It can be expected that a superconducting body, similarly to the one described at Wampler et al. [1], would contain multiple smaller, disjointed domains with complex shapes. The elongation of the superconductor interface in the direction of a magnetic field gradient can lead to a larger effective F m . Additionally, the presence of multiple domains can lead to more complex interactions, such as the appearance of a time-varying angular momentum, and cause the superconductor to rotate in resonance with the changes in magnetic flux. Such effects can be further strengthened by non-uniformities of the external magnetic flux density.
The effects of superconductivity can become non-negligible at larger time scales when considering a superconductor orbiting a larger body. Since the force is usually directed away from the source of magnetic flux, it can be a mechanism counteracting drag and radiation pressure. Additionally, a superconducting bit may significantly change a magnetic environment in its vicinity. It could be especially important in affecting the behaviour of space dust and the formation of larger bodies, if they contained ferromagnetic particles. Generally, though, the possible presence of natural superconductivity at small bodies would have a very small effect on their trajectories. The generation of the natural magnetic fluxes is often connected with the movement of conducting masses and an associated emission of heat. It leads to the increase of equilibrium temperature above the critical temperature in the regions where the magnetic fluxes and their gradients are high enough to generate an appreciable force. This does not exclude the possibility of using artificial cooling and artificial superconductors to allow for a larger force generation. In the case of the terrestrial magnetic field, the effect would still be too weak to be practical. On the other hand, Jupiter's strong magnetosphere and lower temperatures around the planet may make the use of superconducting magnetic sails possible.

Conclusion
The natural existence of a superconducting state in space is possible, and it could have some effect on the behaviour of a small body. However, naturally occurring superconductivity fails to produce a strong magnetic force in viable scenarios. In most cases, a magnetic flux generation is connected with a heat emission. It causes an equilibrium temperature in the vicinity of the magnetic flux source to be above a critical temperature of the superconducting materials, which may be produced and ejected to space by natural processes.
A magnetic flux density and its gradient decrease rapidly with the increasing distance from a body. The derived formulas show the dependence of magnetic force on the distance between the bodies as 1/d 7 . That makes the force at the regions with sufficiently low equilibrium temperature extremely small, to the point where it can be neglected. Any significant effect of superconductivity on the trajectory of a body would require long time scales or a cold source of strong magnetic flux.
The obtained results do not rule out the possibility of using an artificial cooling and hightemperature superconductors to obtain higher forces at specific conditions. It is planned to explore this possibility in the future works on the subject. Superconducting magnetic sails can find use as the method to maintain an orbit by counteracting drag forces. Superconductivity can also play some role in the formation of larger bodies and behaviour of space dust by generating small regions with relatively strong gradients of magnetic field.
Funding Open Access funding provided by CERN.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.