Generation of megatesla magnetic fields by intense-laser-driven microtube implosions

A microtube implosion driven by ultraintense laser pulses is used to produce ultrahigh magnetic fields. Due to the laser-produced hot electrons with energies of mega-electron volts, cold ions in the inner wall surface implode towards the central axis. By pre-seeding uniform magnetic fields on the kilotesla order, the Lorenz force induces the Larmor gyromotion of the imploding ions and electrons. Due to the resultant collective motion of relativistic charged particles around the central axis, strong spin current densities of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim$$\end{document}∼ peta-ampere/\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hbox {cm}^{2}$$\end{document}cm2 are produced with a few tens of nm size, generating megatesla-order magnetic fields. The underlying physics and important scaling are revealed by particle simulations and a simple analytical model. The concept holds promise to open new frontiers in many branches of fundamental physics and applications in terms of ultrahigh magnetic fields.


Scientific RepoRtS
| (2020) 10:16653 | https://doi.org/10.1038/s41598-020-73581-4 www.nature.com/scientificreports/ the implosion, the Lorenz force deflects ions and electrons clockwise and anticlockwise, respectively, gaining azimuthal momentum, as depicted in Fig. 1c. The ion trajectories draw circles with Larmor radii ∼ 0.1 − 1 mm for typical laser and target parameters in MTI. In particular, the envelope of the ion paths forms a nanometerscale hole at the center (hereafter called the "Larmor hole"). Since electrons are negatively charged, the resultant direction of the electron current J eφ is anticlockwise, which is the same as that of the ion current J iφ . Then ultraintense spin currents on the order of 10 15 A cm −2 run around the Larmor hole. Consequently, the currents from the ions and electrons work together to generate MT-order magnetic field B c at the center. Compared with other conventional approaches, the most innovative point of the current concept lies in the geometrically unique plasma flow. A cylindrically converging flow composed of relativistic electrons and ions, which are infinitesimally twisted by the pre-seeded magnetic field in opposite directions, can effectively produce ultrahigh spin currents and consequently, ultrahigh magnetic fields. In addition, the current geometry may be better suited for many practical purposes.
For over 50 years, researchers have strived to realize high magnetic fields. Many approaches have been employed, including high explosives 1, 2 , electromagnetic implosions 3,7 , high-power lasers 10,11 , and Z pinches 5,6 . The principal physical mechanism of these works is based on magnetic flux compression (MFC) using hollow cylindrical structures and pre-seeded magnetic fields. The present scheme also uses a similar physical configuration. However, MTI differs from MFC because the ultrahigh magnetic fields in MTI are generated by the spin currents induced by collective Larmor gyromotions.

Results
Two-dimensional particle simulation. To demonstrate the expected behavior of MTI, we perform 2D (x, y) PIC simulations using the open-source fully relativistic code EPOCH 45 . In this first part of EPOCH simulations (v.4.10.17), we employ rather simple and ideal physical conditions to effectively extract the salient features of the underlying MTI physics. First, the simulation uses the periodic boundary conditions for particles and fields, where the hollow cylindrical volume is placed at the middle of the square computational domain. This configuration simulates collective targets with multiple equally spaced microtubes inside a heated material. We www.nature.com/scientificreports/ set 100 particles/cell for carbon ions and 200 particles/cell for electrons. The lengths for the unit cell size and full span of each side of the square domain are 6.25 nm and 10 µm , respectively. Therefore, the whole computational domain size is 1600 × 1600 mesh 2 . The initial inner radius of the microtube is R 0 = 3 µm. Second, since hot-electron average energy E he.av spans the relativistic regime for the parameters of interest, we use the Maxwell-Jüttner (M-J) distribution 46 rather than the Maxwell-Boltzmann (M-B) distribution for the non-relativistic regime. The M-J distribution defines the hot electron population in terms of the Lorenz factor γ as f (γ ) = γ 2 β �K 2 (1/�) exp(− γ � ) , where β = v/c = 1 − 1/γ 2 and � = T e /m e c 2 with T e , m e , and c being the electron temperature, the electron rest mass, and the speed of light, respectively; K 2 is the modified Bessel function of the second kind. The relation between E he.av and T e significantly differs between the two distributions. That is, E he.av = 3 2 T e for the M-B distribution ( E he.av ≪ m e c 2 ) while E he.av ≃ 3T e − m e c 2 for the M-J distribution ( E he.av ≫ m e c 2 ).
It should be noted that on such an ultrashort timescale as femtoseconds, there is insufficient time for electrons to be thermalized 47,48 . In this sense, employing the M-J distribution, which is characterized by a specific temperature, may not be legitimate. However, high-energy-tail electrons, whose population decreases exponentially with energy, predominantly influence the energy transport and thus the dynamics of the overall system 44 . In fact, both the M-J and M-B distributions have such an exponential dependence in their functional forms. For this reason, employing the M-J distribution is an acceptable choice. Actually, simulations have confirmed that the same value of E he.av yields a similar result for the implosion dynamics and the generation of the magnetic field for both energy distributions. Therefore, E he.av rather than T e is employed below as a principal parameter. Note that we later provide another set of simulation results as a proof-of-principle, using more practical conditions that take the laser−matter interactions into account, where the electron population is not approximated by the M-J distribution. Figure 2 shows the temporal evolution of the normalized densities of ions, ñ i = n i /n i0 , and electrons, ñ e = n e /n e0 , under n e0 = Zn i0 for a fully ionized carbon plasma with A = 12 and Z = 6 . The solid and dashed curves indicate the EPOCH results and the model prediction, respectively. The model is described later. Initially, the inside of the tube is empty, and the remaining volume is filled with uniform ions with T i = 10 eV and n i0 = 1 × 10 23 cm −3 and uniform electrons with E he.av = 5 MeV. The pre-seeded magnetic field B 0 = 4 kT is distributed uniformly over the entire computational domain.
After launching the plasma expansion into a vacuum at τ = 0 , the implosion phase is observed for a period, τ 70 fs (Fig. 2). The implosion velocity of the innermost ions remains nearly constant at v i ≃ 6 × 10 9 cm/s before the cavity collapse. Macroscopically, ions and electrons in the imploding plasma layer move together and maintain charge neutrality. The electron sheath thickness at the plasma/vacuum interface is roughly equal to the local electron Debye length De = (T e /4πn e e 2 ) 1/2 ∼ 150 nm, where e denotes the elementary charge, T e ≈ 1.8 MeV ( E he.av = 5 MeV with the M-J distribution), and n e ≈ 6 × 10 21 cm −3 (Fig. 2). www.nature.com/scientificreports/ Upon cavity collapse, the head group of imploding ions passes the target center at the Larmor hole radius r = R H and expands outward. The mean-free-path of ion-ion collisions is roughly given by ℓ ii ∼ T 2 e /4π n i Z 2 e 4 , which amounts to ∼ 4 cm ≫ R 0 under T e = 2 MeV, Z = 6 , and n i = 10 23 cm −3 . Hence, these ions collisionlessly intersect other ions, which are still imploding toward the center. Meanwhile, the central density increases to the same order as its initial value due to the geometrical accumulation effect (Fig. 2, inset). The Larmor hole radius, R H ≈ 20 nm, indicated at the top of the inset corresponds to the analytical prediction at τ = 90 fs (Eq. 5). The Larmor hole is also seen in the simulation result as the one-humped structure, but the simulated one grows more quickly than the model and expands outward. This is attributed to the fact that a highly compressed ion sphere is created at the center and the strong electrostatic field radially pushes the ions outward. Figure 3 shows snapshots taken from the dominant period of the magnetic field generation, τ ≈ 80 − 120 fs: (a) the normalized ion density ñ i = n i /n i0 , (b) the normalized electron density ñ e = n e /n e0 , (c) the azimuthal electron current J eφ , and (d) the magnetic field B z . Comparing Fig. 3a and b provides insight on how the core plasma develops. The one-humped structure forms in the central region and oscillates at the ion-plasma frequency ω pi = (4πn i Z 2 e 2 /m i ) 1/2 to emit compression waves outward at sound speed c s = (ZT e /m i ) 1/2 . Applying the numbers used in Fig. 3 (i.e., n i = 1 × 10 23 cm −3 and T e = 1.8 MeV) yields c s ≃ 9 × 10 8 cm/s and the cycle τ cyc = 2π/ω pi ≃ 9 fs ( ν ≡ τ −1 cyc ≃ 110 THz), which agree well with the simulation result. According to Ampere's law, c∇ × B = 4πJ +Ė , the azimuthal current distribution, J φ = J eφ + J iφ , directly contributes to the magnetic field B c generated at the center. The azimuthal electron current density J eφ dynamically evolves around the center over the distance approximately equal to the local Debye length De ∼ 100 − 150 nm (Fig. 3c). According to Faraday's law, c∇ × E = −Ḃ , when B c reaches its peak, the displacement current ( ∝ ∂E φ /∂t ) becomes substantially small. Then, B c is given as the sum, Model. Here, we describe the ion dynamics in terms of a semi-analytical model and demonstrate that it forms the basis of the whole system. The time origin matters when comparing the model to the simulation results. To avoid confusion, hereafter, we employ the time variable, t, instead of τ used for the simulation. Suppose that a planar plasma is held at rest in a half-infinitely stretched region −∞ < x ≤ 0 for t ≤ 0 , which is composed of uniform cold ions and hot electrons with densities n i and n e , respectively. We postulate that the plasma is charge neutral, i.e., Zn i = n e . In addition, hot electron temperature T e is assumed to be constant both spatially and temporally due to the high conductivity. Once the boundary between the vacuum and plasma is set free at t = 0 , the plasma begins to expand into the vacuum. The ion motion is governed by the following hydrodynamic system describing the mass and momentum conservation as The two physical ingredients, collisionless ions and isothermal electrons, provide insight to harness Grevich's self-similar solution as a useful approximation and to describe the kinetic behavior of the ions. To accomplish this, a geometrical modification needs to be added to the self-similar solution. The resultant system behaves such that individual fluid elements can penetrate each other in cylindrically converging and diverging processes.
Suppose that an ion with mass m i and ionization state Z is moving at a constant speed v i on the xy-plane in a uniform magnetic field B 0 , which is parallel to the z-axis. The ion draws a circular orbit with a Larmor radius R L = m i v i c/ZeB 0 , where B 0 = |B 0 | . If the position and velocity of the ion are specified at t = 0 to be (x, y) = (R 0 , 0) and (ẋ,ẏ) = (−v i , 0) , respectively, the ion moves on the circle, Here, it is useful to introduce cylindrical coordinates, r = x 2 + y 2 and φ = tan −1 (y/x) . The ion path s along its Larmor circle is measured with the distance from the initial point (x, y) = (R 0 , 0) or with the polar angle θ = sin −1 [(R 0 − x)/R L ] pivoting around the guiding center (x, y) = (R 0 , R L ) , as illustrated in Fig. 4a. It should be noted that Fig. 4a is not to scale. The dimensionless coordinate of the self-similar solution is then redefined in terms of s and θ as ξ = s/c s t = R L θ/c s t . Note that v i and R L are functions of ξ . Consequently, Grevich's selfsimilar solution for a planar system is reformed for a cylindrical system as where v ir and v iφ denote the radial and azimuthal component of the ion velocity, respectively. The reformed set, Eqs. (3) and (4), rigidly satisfies the mass conservation law for a cylindrical geometry such that the factor R 0 /r in Eq.  www.nature.com/scientificreports/ speed of the ion front, i.e., ẋ f → ∞ as ξ → ∞ . Assuming that the plasma expands adiabatically, a reasonable approximation for the dimensionless coordinate of the ion front ξ f is obtained such that the plasma front expands at the speed ẋ f = 2(γ − 1) −1 c s 50 , where γ denotes the adiabatic index. Meanwhile, the self-similar solution based on the isothermal assumption gives the speed of a fluid element at ξ = ξ f as ẋ f = (ξ f + 1)c s . Equating the two speeds gives ξ f = (3 − γ )/(γ − 1) . In particular, the adiabatic index for the relativistic electrons is γ = 4/3 50 , yielding ξ f = 5 . This result well explains ξ f ≃ 5.5 obtained from the simulation (Figs. 2 and 4b). Figure 4b shows the r − t diagram obtained from the reformed self-similar system, under the conditions in Figs. 2 and 3, where the curves correspond to different values of ξ . The time origin of the model corresponding to τ = 12 fs is chosen such that the cavity-collapse timings coincide with each other on the horizontal axes. This can be confirmed by the red dashed curve, which shows the trajectory of the innermost ions at an early stage of implosion according to the EPOCH simulation (Fig. 2). Combining the curves in Fig. 4b and the density profile given by Eq. (3) leads to the dashed curves in Fig. 2 as the model predictions.
The Larmor hole radius R H is obtained from the geometrical consideration under R H ≪ R 0 ≪ R L to be R H ≃ R 2 0 /2R L , which is explicitly rewritten as where ξ c (t) = R 0 /c s t corresponds to the ions passing by the target center at time t. The azimuthal ion current, J iφ = Zen i v iφ , is then given with the help of Eqs. (3) and (4) by The spatial profile of the ion current has a maximum J H (t) at the Larmor hole rim r = R H to spatially decay at the rate r −2 for r ≥ R H . In Eq. (7), the factor R 0 /R H explains the cylindrical accumulation effect. Consequently, the ion current contribution to the central magnetic field follows B ci (t) = (4π/c) With time t, the numerical factor, (ξ c + 1)e −ξ c −1 , in Eq. (7) monotonically increases, and ξ c (t) decreases from its initial value ξ c (τ = 70 fs) = ξ f = 5.5 . Although the factor asymptotically approaches its maximum, e −1 = 0.37 , with t → ∞ or ξ c → 0 , a cut-off value of ξ c exists for the physical reason described below.
After the duration �τ , the core periodically emits outgoing density waves at a frequency ω pi . These waves carry a portion of the central plasma energy, as seen in the double-humped structure of the ion density profile for τ = 110 − 115 fs in Fig. 3a. In fact, B c (t) begins to decay coherently with the first emission of the density wave (Fig. 3a,d). Since the model does not consider this emission process, we estimate the maximum magnetic field by limiting the growth at the peak time τ peak = τ c + �τ . This corresponds to ξ c ∼ 3 . Recalling the observed constancy, B c ≃ 4B ci , leads to its maximum value B c.max as Thus, B c.max is proportional to the total ion flux emitted from the inner surface of the microtube, i.e., B c.max ∝ Zn i0 c s R 0 (recall c s ∝ √ ZT e /A). Figure 5 shows the results of about three dozen EPOCH simulations (solid circles) for B c.max as a function of defined in Eq. (8). The straight black line denotes the model prediction. Each simulation result corresponds to a subset with the key parameters, (B 0 , n i0 , R 0 , E he.av ) . Their composition is chosen rather randomly over the ranges, B 0 = 1 − 10 kT, n i0 = 5 × 10 22 − 2 × 10 23 cm −3 , R 0 = 1 − 3 µm , and E he.av = 5 − 15 MeV, while A = 12 and Z = 6 are fixed. Despite the random choices, the overall simulation results are smoothly linked in a systematic manner by the dashed curves parameterized by B 0 . This demonstrates the physical significance of the parameter as an essential measure of the magnetic-field generation in the present scheme. There is a threshold relation between B 0 and such that the simulation results and the model line agree well. For example, for B 0 6 kT and 2.5 , the model well reproduces the overall behavior of the simulation results. The physical reason why these curves overlap with the model line is as follows. Although the individual trajectories of charged particles are deformed to radially shift outward due to higher B 0 and smaller Larmor radius R L (Fig. 4a), the integrated currents and consequently the magnetic field, J iφ dr ∝ J eφ dr ∝ B c , are unchanged.
With B 0 = 4 kT, the difference between the simulation and the model begins to increase for 1.6 . Moreover, for B 0 < 3 kT, the behaviors of the simulation results are unpredictable by the model. In particular, with B 0 = 2 kT, the temporal evolution of the system becomes unstable. Although it initially behaves in the reverse polarity regime ( B c.max < 0 ), it eventually turns into the forward polarity regime ( B c.max > 0 ). This phenomenon is labeled as "polarity switching" in Fig. 5. The polarity switching suddenly occurs on the timescale of several femtoseconds, when the electron current distribution surrounding the target center evolves very quickly in a complex manner. One of the potential causes for this phenomenon seems to be the existence of an electron-rich www.nature.com/scientificreports/ space at the center, which is omitted in the model assuming that R H is so large (Eq. 5) that electrons are evacuated from the Larmor hole.
Practical simulation with laser-plasma interaction. In this section, we perform proof-of-principle EPOCH simulations to demonstrate that strong magnetic fields can still be produced, when the uniform hot electron population previously assumed by the M-J distribution is replaced with a realistic laser-plasma interaction 45 . These simulations reproduce salient features of the MTI process described in the previous sections. Note that although the EPOCH code used here is somewhat customized with the same core modules as recent releases, the results provided below should be reproducible using the latest release. We consider an isolated target, which consists of an initially cold, fully ionized charge-neutral carbon-electron plasma with an initial ion density n i0 = 3 × 10 22 cm −3 and electron density n e0 = 6n i0 . We have reduced the target density by a factor of ∼ 3 relative to the case given in Figs. 2 and 3 to mitigate the computational cost associated with PIC simulations of the laser−plasma interaction. The target's outer cross-section is a square with 12-µm-long sides. The inner radius of the microtube is R 0 = 3 µm . This target is irradiated on each of the four outer sides by a large-spot (spatially plane wave) laser pulse with a wavelength of L = 0.8 µm , a total duration ( sin 2 temporal shape in |E|) of τ L = 100 fs, and a peak intensity of I L = 10 21 W/cm 2 . The pulses are co-timed such that the peaks of all four pulses interact with the target surface simultaneously. The plasma is modeled with a resolution of 100 cells/µm and 200 particles/cell for carbon ions and 400 particles/cell for electrons, using the cubic B-spline particle shape included in EPOCH. The full size of the simulation box is 22 µm × 22 µm.
In practical laser-driven configurations, the laser-plasma interaction causes the hot electron population driving the implosion to depart from the conditions of spatial uniformity and temperature isotropy, which are assumed in the derivation of the maximum magnetic field given in Eq. (8). The energetic electron spectrum has multiple energy components and differs from a single-temperature M-J distribution (Fig. 6a). Despite this departure from the conditions assumed in the previous sections, we still observe strong magnetic field generation with similar features to the single-temperature case. Due to the spatial non-uniformities, strong magnetic fields with spatially varying signs can be generated in small regions of the implosion volume even without the presence of any seed magnetic field. However, the addition of a seed magnetic field ( B 0 = 6 kT) increases both the maximum amplitude of the magnetic field produced and the spatial volume over which it maintains the same sign (Fig. 6b). Figure 7 shows the detailed temporal evolution of the physical quantities, which are summarized in Fig. 6b. The magnetic field generation in the laser-driven case occurs in two stages, where the majority of the magnetic field generation occurs during the first stage. During the implosion, substantial anisotropy in the ion current crossing through the center of the microtube generates a strong magnetic field around the center ( r < 0.3 µm , www.nature.com/scientificreports/ Stage 1 in Fig. 7a,c). In the simulations, the magnitude of this magnetic field increases by a factor of 2 upon applying B 0 = 6 kT through the MTI process (Fig. 7a). Later, this central magnetic field is further amplified by electrons undergoing E × B-directed motion as the central ion population explodes outward (Stage 2 in Fig. 7a,c). Unlike Stage 1, which occurs over approximately 30 fs and agrees well with the MTI process described earlier, Stage 2 can persist for well over 100 fs. Figure 7 does not capture the end of this stage due to the computational cost of these simulations. However, simulations performed with plastic targets suggest the magnetic field can be  Detailed temporal evolution of the magnetic field B z and the normalized ion density n i /n i0 of laserdriven MTI. Summary plots are given in Fig. 6. (a,b) results with B 0 = 6 kT. www.nature.com/scientificreports/ slowly amplified over hundreds of femtoseconds. In the presence of the 6-kT seed, the first stage, which includes the MTI amplification process, generates a ∼ 100 kT magnetic field over r 0.3 µm in approximately 20 fs, while the second stage amplifies this magnetic field to ∼ 120 kT over approximately 50 fs and increases the size of the central spot to a radius of ∼ 0.5 µm.

Discussion
We here briefly discuss laser requirements for MTI. To achieve MT-order magnetic fields experimentally, a rough estimate assuming a pulse duration of ∼ 30 fs suggests that a laser system with a pulse energy of 0.1 − 1 kJ and a total power of 10 − 100 PW is required. Such high-power laser performance is accessible by today's laser technology 51,52 . Meanwhile, fundamental studies, including proof-of-principle experiments, should be feasible using substantially smaller laser systems. Unlike for ultrathin targets with nm-scale thicknesses, MTI targets are significantly less sensitive to laser contrast due to the micron-thick wall, while beam co-timing should be within 10 − 20 fs for implosions with a timescale of ∼ 100 fs. Detecting MT-order magnetic fields inside a plasma presents a challenge for conventional techniques that rely on charged particle sources. In anticipation of achieving ultrahigh magnetic fields, there have been efforts to develop other techniques to infer the existence of strong B-fields inside a dense plasma in use of, for example, an XFEL photon beam with Faraday rotation effect 17,20,24 and spin-polarized neutrons 53 .
We roughly estimate the minimum number of beams n B from a uniformity point of view. For simplicity our discussion here is limited to the cross-sectional dimensions. Nonuniformity of the imploding ion front is directly influenced by that of hot electron density, which comprises the local electron sheath (Fig. 1b). Hot electrons produced on the laser-irradiated surface go back and forth between the target surface and the cavity wall with an in-between distance R = R 1 − R 0 , where R 1 is the initial target radius. This hot electron transport is regarded as a kind of random-walk diffusion process. The nonuniformity on the outer surface should diminish on the cavity wall via this diffusion process along the lateral direction. As is well documented, the diffusion distance is given by L d ∼ √ N R , where N ∼ c�t/�R stands for the reflection number between the two surfaces during the implosion time t . Consequently, n B 2πR 0 /2L d ∼ 2 (i.e., two-sided illumination) when employing for example R 0 = 3 µm , �R = 2 µm , and t ∼ 50 fs (Fig. 2). It should be noted that controlling both the temporal and spacial profiles of the incident laser pulses also plays a crucial role to improve the implosion uniformity 54,55 .
In summary, we propose a novel concept called MTI, which produces MT-order magnetic fields using intense laser pulses. Key physical elements of MTI are imploding ion fluxes with quasi-relativistic speeds and the resultant ultrahigh spin currents running around the nanometer-scale Larmor hole at the center. The spin currents are due to the collective motion of the imploding ions and the accompanying relativistic electrons. The pre-seeded magnetic field B 0 significantly influences the magnetism of the plasma. For example, the forward (reverse) polarity appears in the domain of higher (lower) B 0 . Polarity switching is an extraordinary phenomenon, which requires further investigation. The scaling law for the maximum magnetic field B c.max is obtained as a function of B 0 and the total ion flux emitted from the inner surface of microtube . With the realistic laser-plasma interaction taken into account, strong magnetic field generation has been demonstrated as a proof-of-principle of MTI.