# Inverse effect of magnetostriction in magnetoelectric laminates

M. Auslender  E. Liverts  B. Zadov  A. Elmalem  A. Zhdanov  A. Grosz  E.
Paperno  Ben-Gurion University of the Negev, P.O. Box 653, Beer-Sheva 84105,
Israel

###### Abstract

We introduce the notion of inverse effect of magnetostriction for
magnetostrictive-piezoelectric heterostructures and study this effect
theoretically and experimentally. It is shown that the inverse effect of
magnetostriction may crucially contribute to the mechanism of magnetoelectric
coupling. It is shown that the studied effect essentially modifies the
saturation magnetostriction of the whole structure as compared to its magnetic
phase bulk and also induces an additional magnetic anisotropy. Our
consideration provides useful insight into the fundamental issue of strain-
mediated magnetoelectric coupling. Understanding this effect may lead to its
utilization in original experimental concepts and the improvement of the ME
coupling.

###### pacs:

Valid PACS appear here

Magnetoelectric (ME) laminates that comprise bonded parallel magnetostrictive
and piezoelectric layers are important subfamily of two-phase multiferroics W.
Eerenstein, N. D. Mathur, J. F. Scott . In these heterostructures, the
indirect strain-mediated ME coupling results from an interplay of
magnetostriction in the magnetic layer, piezoelectricity in the piezoelectric
layer, and the bonding. So far, inverse effect of magnetostriction was meant
as a change in the magnetic state of magnetostrictive media subjected to
external stresses, being strong when the magnetoelastic and magnetic
anisotropy energies are comparable bozorth ; Chikazumi ; magnet ; Miyauchi
(1987); Livingston (1982); Zheng . Recently, similar effect of residual
stresses in a magnetostrictive film on a compliant substrate was observed Li .
In the ME laminates subjected to an external magnetic field, internal stresses
emerge as a result of mechanical compliance of the phases. As regards the
magnetic phase, these stresses will drive the inverse effect of
magnetostriction in the usual sense. In addition due to the hidden stresses,
the whole laminate will deform upon increasing the magnetic field, quite
differently than free standalone magnetic layer. We define these two phenomena
altogether, as the inverse effect of magnetostriction for ME laminates. Though
the strain-mediated ME coupling in the laminates has been intensively studied,
W. Eerenstein, N. D. Mathur, J. F. Scott ; srin ; Nan, M. I. Bichurin, S.
Dong, D. Viehland, and G. Srinivasan (2008) and references therein, no
specific attention has been paid to the laminate derived inverse
magnetostriction effect. Even more, we believe that just missing the inverse
effect of magnetostriction from the scope is a physical reason behind drastic
theoretical overestimation of the ME coupling as compared to its experimental
value Laletin et.al (2005); liverts1 . The present work reports the study of
the inverse effect of magnetostriction for the ME laminates. Restricting
ourselves to a range of magnetic field near the magnetic saturation allows us
to combine physical rigor with relatively simple treatment. Our consideration
concerns the polycrystalline magnetic phase, which is consistently taken into
account, though the effect under study works also for crystalline magnetic
one. We prove that when the phase has low magnetic anisotropy, the inverse
effect of magnetostriction is crucial and may change the ME coupling by an
order of magnitude as compared to the estimate obtained from piezomagnetic
coefficient for the magnetic phase alone. The proposed account of the inverse
effect of magnetostriction also provides emergent experimental concepts in the
field of ME heterostructures.

Magnetostrictive polycrystalline layers, whose structural components are
sufficiently small single-crystal grains, may be regarded as an isotropic body
in terms of both elastic and magnetostrictive properties, especially in the
range close to the magnetization saturation (since we are dealing with the
strains and magnetization distribution over the spatial regions that are large
compared to the grains size) landau_el ; Chikazumi ; magnet . In the cases
where the anisotropy of the elastic properties is relatively small, and
magnetostriction constants are of the same sign (e.g. Ni), the
magnetostriction constants and elastic modules of the isotropic body can be
accurately found from those of the single crystal by straightforward
calculating the strains in each grain and then averaging them over the grains
landau_el ; magnet .

For a strained cubic crystal, the magnetoelastic energy expressed for each
domain with a local frame $x,\;y,\;z$ in terms of the strain $u_{ik}$,
magnetization $\mathbf{M}$, $\mathbf{m}=\mathbf{M}/M_{s}$ ($M_{s}$ is the
saturation magnetization), and two independent coefficients $a_{1}$, $a_{2}$,
is given by

$\mathcal{U}_{\rm
m-el}=a_{1}\left(u_{xx}m_{x}^{2}+u_{yy}m_{y}^{2}+u_{zz}m_{z}^{2}\right)+a_{2}\left(u_{yz}m_{y}m_{z}+u_{xz}m_{x}m_{z}\\\
+u_{xy}m_{x}m_{y}\right).$ (1)

The change in the deformation resulting from a change in the magnetization
direction and from persistent stress $\sigma_{ik}$ are obtained by minimizing
energy $\mathcal{U}^{m}_{\rm el}+\mathcal{U}_{\rm m-el}$, where
$\mathcal{U}^{m}_{\rm el}$ is the elastic energy involving three independent
elastic constants ($s^{m}_{11}$, $s^{m}_{12}$, $s^{m}_{44}$), is equal to

$u_{xx}=s^{m}_{11}\sigma_{xx}+s^{m}_{12}(\sigma_{yy}+\sigma_{zz})+\frac{3}{2}\lambda_{100}\left(m_{x}^{2}-\frac{1}{3}\right),\;u_{xy}=\frac{s^{m}_{44}}{2}\sigma_{xy}+3\lambda_{111}m_{x}m_{y},$
(2)

where $s^{m}_{44}a_{2}=-3\lambda_{111}$,
$2(s^{m}_{12}-s^{m}_{11})a_{1}=3\lambda_{100}$, and similarly for the other
components. According to landau  the above strain tensor presents all the
magnetostriction effects in the cubic crystals.

Let $n_{x}$, $n_{y}$, $n_{z}$ be the cosines of a fixed direction $\mathbf{n}$
in the local cubic crystallographic frame of a crystallite, which serves as
the axis $1$ of a laboratory frame with the axes $1$, $2$, $3$. Then, we adopt
the following transformation between the local and laboratory frame

$x_{1}=xn_{x}+yn_{y}+zn_{z},\;\;x_{2}=x\frac{-n_{y}\cos\psi+n_{x}n_{z}\sin\psi}{\sqrt{1-n_{z}^{2}}}+y\frac{n_{x}\cos\psi+n_{y}n_{z}\sin\psi}{\sqrt{1-n_{z}^{2}}}-z\sqrt{1-n_{z}^{2}}\sin\psi.$
(3)

The expression for $x_{3}$ can be obtained from the equation for $x_{2}$ under
the replacement $\psi\rightarrow\psi-\pi/2$. Assume that within the
magnetostrictive layer the strains, stresses, and magnetization are nearly
uniform and close to their saturation values. Then using Eq. (3), averaging
the strains in Eq. (2) over all the directions of $\mathbf{n}$ (using the
following computations $\langle n_{i}^{4}\rangle=1/7$, $\langle
n_{i}^{8}\rangle=1/9$ etc.), $0\leq\psi\leq 2\pi$, and comparing the results
with the appropriate expressions for an isotropic medium, we find that the
elastic properties can be described by the Young’s modulus $E$ and Poisson’s
ratio $\nu$, given by

$\frac{1}{E}=\frac{3}{5}s^{m}_{11}+\frac{1}{5}\left(2s^{m}_{12}+s^{m}_{44}\right),\quad
2\nu=-\frac{2s^{m}_{11}+8s^{m}_{12}-s^{m}_{44}}{3s^{m}_{11}+2s^{m}_{12}+s^{m}_{44}},$
(4)

and the magnetostriction can be described by a single coefficient
$\lambda_{s}$, as in truly isotropic material.

To describe the inverse effect of magnetostriction in a bonded thin
magnetostrictive-piezoelectric layered structure, we will define the in-plane
coordinate axes as $1,\;2$ and the off-plane axis as $3$. Under an external
magnetic field along the $1$ direction ($H_{1}=H,\;H_{2}=H_{3}=0$), the
magnetostrictive layer shrinks or expands along the $1$, $2$ directions. If
there is no slippage or fracture in the bonded trilayer, the magnetostrictive
strain applies an in-plane constraint on the compliant piezoelectric layer
which also can either shrink or stretch to minimize its energy
$\mathcal{U}^{p}_{\rm el}+\mathcal{U}_{\rm e-el}$ and can help the
piezoelectric to have induced polarization along the 3 axis according to its
phase diagram. Here $\mathcal{U}_{\rm e-el}$ is the usual electro-elastic
energy of the piezoelectric, which we take in a strain-charge form landau ;
mason . For the 1-2 plane considered to be infinite, the uniform stresses
corresponding to the uniform strain, are then imposed on each of the two
phases to maintain displacement compatibility and equilibrium. These stresses
are expressed via the layer’s thicknesses and piezostrictive,
magnetostrictive, and elastic properties; the off-plane stress components are
neglected due to the small thickness and the fact that the trilayer structure
prevents bending. Taking into account that the forces per unit area exerted by
the stresses in the magnetostrictive layers are balanced by those caused by
the stresses in the piezoelectric and considering open-circuit condition, we
obtain at the magnetic saturation

$\sigma_{11}=-\frac{t_{p}\lambda_{s}}{2(s_{11}-s_{12})}\frac{2s_{11}+s_{12}+3t_{m}d_{31}^{2}/\varkappa}{s_{11}+s_{12}-2t_{m}d_{31}^{2}/\varkappa},\;\;\;\sigma_{22}=\frac{t_{p}\lambda_{s}}{2(s_{11}-s_{12})}\frac{s_{11}+2s_{12}+3t_{m}d_{31}^{2}/\varkappa}{s_{11}+s_{12}-2t_{m}d_{31}^{2}/\varkappa},$
(5)

where $t_{p}$ is the piezoelectric thickness, $t_{m}=1-t_{p}$,
$s_{11}=t_{p}E^{-1}+t_{m}s^{p}_{11}$, $s_{12}=-\nu
t_{p}E^{-1}+t_{m}s^{p}_{12}$, $s^{p}_{ik}$, $d_{31}$ and $\varkappa$ are the
piezoelectric elastic constants, piezoelectric coefficient, and out-of-plane
dielectric permittivity, respectively. As it recently has been argued liverts
and also confirmed numerically zadov  by a finite element method, the
expressions in Eq.(5) are leading asymptotic in terms of the lateral-to-
transverse dimensions ratios tending to infinity. The simulations zadov  for
finite-lateral sizes’ structures have shown deviation from Eq.(5) only in
transition layers near the edges, where the physical fields are strongly
inhomogeneous.

One of the outcomes, most important for applications, is that a significant
change in the averaged deformation occurs in the laminate as compared to the
deformation in the bare magnetostrictive phase. The average strains at $H\gg
M_{s}$, resulting from the magnetostriction, are $\langle
u_{11}\rangle=(\sigma_{11}-\nu\sigma_{22})E^{-1}+\lambda_{s}$ and $\langle
u_{22}\rangle=(\sigma_{22}-\nu\sigma_{11})E^{-1}-\lambda_{s}/2$. Then the
calculated magnetostriction $\lambda_{\theta}$ at an angle $\theta$ to the
magnetization is given by

$\lambda_{\theta}=\frac{\delta
l_{\theta}}{l}=\frac{1-\nu}{2E}(\sigma_{11}+\sigma_{22})+\frac{\lambda_{s}}{4}+\left[\frac{1+\nu}{2E}(\sigma_{11}-\sigma_{22})+\frac{3\lambda_{s}}{4}\right]\cos
2\theta.$ (6)

The measurement of $\lambda_{\theta}$ vs $\theta$ allows one to extract the
lamination stresses and compare them with those theoretically predicted by
Eq.(5).

At elastic equilibrium between the bonded magnetostrictive and piezoelectric
layers, the total energy of the former
$\mathcal{U}^{m}\left(\mathbf{m},\mathbf{H}\right)$ is eventually obtained by
substituting the equilibrium strains, see Eq.(2), into the expression of
$\mathcal{U}^{m}_{\rm el}+\mathcal{U}_{\rm m-el}$. This yields
$\mathcal{U}^{m}=\mathcal{U}^{m}_{0}+\Delta\mathcal{U}^{m}$, where
$\mathcal{U}^{m}_{0}$ is sum of the Zeeman energy and the energy of
crystallographic magnetic anisotropy of the bare magnetostrictive phase at
zero strain, and $\Delta\mathcal{U}^{m}$ is given by

$\Delta\mathcal{U}^{m}\left(\mathbf{m}\right)=\Delta\mathcal{U}^{m}_{c}\left(\mathbf{m}\right)-\frac{3}{2}\lambda_{100}(m_{x}^{2}\sigma_{xx}+m_{y}^{2}\sigma_{yy}+m_{z}^{2}\sigma_{zz})-3\lambda_{111}(m_{x}m_{y}\sigma_{xy}+m_{x}m_{z}\sigma_{xz}+m_{y}m_{z}\sigma_{yz}).$
(7)

Here $\Delta\mathcal{U}^{m}_{c}$ is the cubic anisotropy energy induced by the
magnetostriction in a bulk ferromagnetic at zero stress (coinciding with
expressions given in the literature, see e.g. Chikazumi ), while the rest
terms represent the magnetic anisotropy due to the stress induced by the
lamination, which is of our primary interest. When the stress is given and is
independent of $\mathbf{m}$, the stress-induced anisotropy is of uniaxial
type. Although
$|\Delta\mathcal{U}^{m}\left(\mathbf{m}\right)|<|\mathcal{U}^{m}_{0}\left(\mathbf{m},0\right)|$,
the induced anisotropy is not negligible. In what follows, we consider the
magnetization curve of the laminate near the magnetic saturation. To this end,
we search the minimum of the above magnetic energy (normalized to
$\mu_{0}M_{s}^{2}$), keeping strict condition $|\mathbf{m}|=1$ with the use a
perturbation theory in small parameter $\varepsilon=M_{s}/H$. Let us write the
minimized function in the form

$\frac{\beta}{4}\left(m_{x}^{4}+m_{y}^{4}+m_{z}^{4}\right)+\frac{\Delta\mathcal{U}^{m}\left(\mathbf{m}\right)-\Delta\mathcal{U}^{m}_{c}\left(\mathbf{m}\right)}{\mu_{0}M_{s}^{2}}-\frac{1}{\varepsilon}\mathbf{n}\cdot\mathbf{m},$
(8)

where $\beta\mu_{0}M_{s}^{2}=2(-K_{1}+\Delta K_{1})$, $K_{1}$, and $\Delta
K_{1}$ are the cubic anisotropy constants at zero strain and the
magnetostriction correction to it, respectively. To obtain the equations
determining $\mathbf{m}$, the conditions for the above minimum is used, which
gives

$m_{x}=\frac{n_{x}}{\lambda}-\varepsilon\left[\beta
n_{x}^{3}-3\frac{\lambda_{100}n_{x}\sigma_{xx}+\lambda_{111}(n_{y}\sigma_{yy}+n_{z}\sigma_{zz})}{\mu_{0}M_{s}^{2}}\right]+O(\varepsilon^{2}),$
(9)

here $\lambda=1+O(\varepsilon)$ is the Lagrange multiplier and similar
equations hold for the components $m_{y,z}$. The measurable quantity is the
magnetization projection $M_{1}=\mathbf{M}\cdot{\mathbf{n}}$ averaged over the
all possible orientations of the individual grains. Averaging Eq.(9) yields

$\frac{\langle
M_{1}\rangle}{M_{s}}=\left\langle\sqrt{1-|\mathbf{m}\times{\mathbf{n}}|^{2}}\right\rangle=1-\frac{1}{2}\left\langle|\mathbf{m}\times{\mathbf{n}}|^{2}\right\rangle+O(\varepsilon^{3}).$
(10)

To accomplish our task, we express the saturated in-crystallite (local) stress
components via their non-zero components in the laboratory frame,
$\sigma_{11}$ and $\sigma_{22}$, using the transformation presented by Eq.
(3), substitute the results into Eq. (10) and average it over all the
directions of $\mathbf{n}$ and the angle $\psi$. This procedure, up to leading
terms in the magnetostriction parameters, results in the law of approach to
the magnetic saturation

$\langle M\rangle\simeq
M_{s}\left[1-\frac{1}{105}\frac{H_{A}^{2}(\sigma)}{H^{2}}\right],\;H_{A}^{2}(\sigma)=2M_{s}^{2}\left[\beta^{2}-\frac{3\beta}{\mu_{0}M_{s}^{2}}(2\sigma_{11}-\sigma_{22})\Delta\lambda\right],\;\Delta\lambda=\lambda_{100}-\lambda_{111}.$
(11)

To verify the modeling of the inverse magnetostriction effect developed above,
we bonded several laminates. Our raw materials were $99.95\%$ purity
polycrystal Ni and APC-844 ceramic lead zinc titanate (PZT), obtained from
Alfa Aesar and APC, respectively. Two $20\times 20$ mm2 area samples with
$0.5$ thick Ni and $0.83$ mm thick laminates were used for magnetostriction
measurements, and another pair of $5\times 2$ mm2 area samples with $0.25$
thick Ni and $0.58$ mm thick laminates, for magnetization measurements. The
magnetostriction measurements were carried out using two-axis SR-4 strain
gauges manufactured by Vishay Micro-Measurements. The magnetization curve at
room temperature were measured on a PPMS machine. Fig. 1 shows the saturated
magnetostriction $\lambda_{\theta}$ and $\lambda_{\pi/2-\theta}$ vs $\theta$.
With material parameters $E_{11}=76~{}\rm{GPa}$, $d_{31}=-109~{}\rm{pC/N}$,
$\varkappa=13.3~{}\rm{nF/m}$ (PZT), $E=160~{}\rm{GPa}$, $\nu=0.35$, and Eq.
(5) we calculated the lamination stresses and used them for obtaining the
theoretical curves. It is seen that a linear function of $\cos 2\theta$ as
predicted by Eq.(6) fits good the measured dependences with the extracted
$\lambda_{s}=-3.4\times 10^{-5}$.

Figure 1:  Magnetostriction as a function of angle $\theta$ for Ni and a
Ni/PZT/Ni laminate. The curves are theoretical predictions for: solid –
$\lambda_{\theta}$, dashed - $\lambda_{\pi/2-\theta}$; $1$, $2$ – the
Ni/PZT/Ni laminate, and $3$, $4$ – Ni. The markers represent the experimental
data: squares – Ni/PZT/Ni laminate, circles – Ni.

Fig. 2 shows a comparison of the calculated $\langle M\rangle$ vs $H$ curves
in high-$H$ range against the experimental $M(H)$ data. For the laminate
measured, we plugged into Eq. (11) the stresses $\sigma_{11}=35~{}\rm{MPa}$
and $\sigma_{22}=-16~{}\rm{MPa}$ estimated using the numerical data zadov .
For the fitting parameters $\beta=0.14$, $M_{s}=458~{}\rm{kA/m}$ and
$\Delta\lambda=-4.6\times 10^{-5}$ the theory agrees well with the experiment,
as clearly seen from Fig.2.

Figure 2:  The magnetization curves in a high-$H$ range. The curves represent
theoretical dependences and the markers represent the experimental data. The
solid curve and squares correspond to the Ni/PZT laminate and the dashed curve
and circles correspond to Ni alone.

We have considered the inverse effect of magnetostriction in the ME laminates
containing a magnetostrictive cubic ferromagnetic polycrystal. The lamination
develops the stresses in both the magnetostrictive and piezoelectric layer,
and the stresses lead to the effective uniaxial magnetic anisotropy, see Eq.
(7). It is worth emphasizing that, similarly to ferromagnetic polycrystal
landau ; Chikazumi , $\langle M\rangle$ of the laminates approaches the
magnetic saturation in the inverse power-two law manner, while the lamination
changes the value of the effective anisotropy field $H_{A}(\sigma)$, see Eq.
(11). On the contrary, an initial susceptibility-adjusted Langevin curve Zheng
; Li  for $\langle M\rangle$ vs $H$ predicts unrealistically fast exponential
saturation. For a Ni/PZT laminate, we have shown that, as compared to a bulk
ferromagnetic, the magnetostriction of the laminate decreases, see Eq. (6) and
Fig. 1, and the magnetic anisotropy field increases, see Eq. (11) and Fig. 2.

The considered effect provides helpful insight into and adds additional
feature to the mechanism of the strain-mediated ME coupling. Indeed in our
notations, the thumb rule estimation of the ME coupling contains the factor
$\left(\lambda_{\theta}-\lambda_{\pi/2-\theta}\right)/H_{A}$. For Ni as the
magnetostrictive phase, the saturated magnetostriction and magnetic anisotropy
field decreases and increases, respectively, in a concurrent manner. This
explains why the measured ME coupling is smaller than the theoretical
estimation of it obtained without including the inverse effect of
magnetostriction, as noted in Introduction. Also in the present framework, the
observed shift of the magnetic field of the optimal ME response to higher
values relative to one estimated from the magnetization curve of the bulk
ferromagnetic phase finds natural explanation. The outcomes of our theory and
their experimental tests are reported here for a range near the magnetic
saturation, but the opposite range of magnetic fields is treated similarly.
Our study leads us to conclude that tailoring the stresses may essentially
improve the strain-mediated ME coupling by compensating the inverse effect of
magnetostriction. The account of this effect may also help for engineering the
ME heterostructures with low magnetic bias and reduced adverse effects of
domain structure, e.g. Barkhausen jumps.

The authors are grateful to I. Felner and M. Tsindlekht from Hebrew University
for their help in the magnetization curves measurements.

## References

  * (1) W. Eerenstein, N. D. Mathur, J. F. Scott, Nature, 442, Issue: 7104, pp. 759-765, (2006)
  * (2) R. M. Bozorth, Ferromagnetism, p. 847, (Princeton, N.J. : D. Van Nostrand, New York, 1951).
  * (3) S. Chikazumi, “Physics of Ferromagnetism,” (New York: Oxford University Press Inc., 1997).
  * (4) Magnetism: Fundamentals, eds. E. Du Tremolet De Lacheisserie, D. Gignoux, M. Schlenker, (New York : Springer, 2005).
  * Miyauchi (1987) T. Miyauchi, T. Kira, M. Yoshikawa, _IEEE T. Magn._ , 23, (5), pp. 2497-2499, (1987)
  * Livingston (1982) J.D. Livingston, _Phys. Status Solidi A: Applied Research_ , 70, 2, pp. 591-596 (1982)
  * (7) X. J. Zheng, X. E. Liu, J. Appl. Phys, 97, 053901, (2005)
  * (8) X.Y. Lu, and H. Li, _J. Magn. Magn. Mater._ , 322, (15), pp. 2113-2116, (2010)
  * Nan, M. I. Bichurin, S. Dong, D. Viehland, and G. Srinivasan (2008) C.-W. Nan, M. I. Bichurin, S. Dong, D. Viehland, and G. Srinivasan, _J. Appl. Phys_ , 103, 031101 (2008)
  * (10) G. Srinivasan, E. T. Rasmussen, B. J. Levin, and R. Hayes, Phys. Rev. B 65, 134402 (2002)
  * Laletin et.al (2005) V.M. Laletin, N. Paddubnaya, G. Srinivasan, C.P. De Vreugd, M.I. Bichurin, V.M. Petrov and D.A. Filippov _Appl. Phys. Lett._ 87, 222507 (2005).
  * (12) E. Liverts A. Grosz, B. Zadov, M. I. Bichurin, Y. J. Pukinskiy, S. Priya, D. Viehland, E. Paperno, _J. Appl. Phys._ , 109, 07D703, (2011).
  * (13) L. D. Landau and E. M. Lifshitz, Theory of elasticity (Pergamon, Oxford, 1970)
  * (14) L. D. Landau and E. M. Lifshitz, Electrodynamics of Continuous Media (Pergamon, Oxford, 1960)
  * (15) D. A. Berlincourt, D. R. Curran, H. Jaffe, in _Physical Acoustics_ , ed. W. P. Mason, Vol. I, part A, p.190, (Academic, New York, 1964)
  * (16) E. Liverts, M. Auslender, A. Grosz, B. Zadov, M. I. Bichurin, and E. Paperno, _J_. Appl. Phys., 107, 09D914, (2010)
  * (17) B. Zadov, A. Elmalem, E. Paperno, I. Gluzman, A. Nudelman, D. Levron, A. Grosz, S. Lineykin, and E. Liverts, _Adv. Condens. Matter Phys._ , 383728 (2012)


# In-Hand Pose Estimation and Pin Inspection for Insertion of Through-Hole
Components

Frederik Hagelskjær and Dirk Kraft  All authors are with The Maersk Mc-Kinney
Moller Institute, University of Southern Denmark.
<EMAIL_ADDRESS>
###### Abstract

The insertion of through-hole components is a difficult task. As the
tolerances of the holes are very small, minor errors in the insertion will
result in failures. These failures can damage components and will require
manual intervention for recovery. Errors can occur both from imprecise object
grasps and bent pins. Therefore, it is important that a system can accurately
determine the object’s position and reject components with bent pins.

By utilizing the constraints inherent in the object grasp a method using
template matching is able to obtain very precise pose estimates. Methods for
pin-checking are also implemented, compared, and a successful method is shown.

The set-up is performed automatically, with two novel contributions. A deep
learning segmentation of the pins is performed and the inspection pose is
found by simulation. From the inspection pose and the segmented pins, the
templates for pose estimation and pin check are then generated.

To train the deep learning method a dataset of segmented through-hole
components is created. The network shows a 97.3 % accuracy on the test set.
The pin-segmentation network is also tested on the insertion CAD models and
successfully segment the pins.

The complete system is tested on three different objects, and experiments show
that the system is able to insert all objects successfully. Both by correcting
in-hand grasp errors and rejecting objects with bent pins.

##  I Introduction

This paper presents a method for post-grasp pose estimation and pre-insertion
inspection. The system performs in-hand pose estimation and pin inspection of
through-hole components. An example of correcting an erroneous grasp is shown
in Fig. 1.

The insertion of objects is an essential task in robotics. For insertions to
succeed, very precise position information of the object is generally
required, and the object should not be damaged. Several methods have been
developed to obtain precise poses of objects. One approach is to focus on
restricting the object [1] or aligning the object with mechanical solutions
[2, 3]. The challenge with such methods is that new hardware is needed when
new objects are introduced, thus limiting the flexibility of the system.

A different approach is visual pose estimation of the object. This is
generally done pre-grasping and has shown very good results in recent years
[4, 5]. As these methods are used pre-grasping, any errors introduced by the
grasping is not compensated. While using these methods post-grasping is
possible, they are challenged by the occlusion presented by the fingers and a
general lack of precision needed for very accurate manipulation [6].

Methods have been developed specifically for post-grasp pose estimation [4,
7], these methods use depth information that is difficult to obtain with small
objects and thus these methods typically focus on larger objects. This paper
introduces a method using 2D data for in-hand pose estimation of very small
objects. The method uses stable poses introduced implicitly by the grasping to
return very precise pose estimations.

For the manipulation to be successful, the objects should also be in
satisfactory condition. If parts are damaged and the object does not
correspond to the rigid CAD model, the manipulation can fail.

The focus of this paper is the insertion of through-hole components. As the
tolerances of the holes are tiny, the precision of the in-hand pose must be
very high, and the pins cannot be bent.

(a) Pose estimation of an object with large angular error compared with the
finger.

(b) During insertion the finger is rotated to compensate for the large angular
error.

Figure 1:  Example of the system’s ability to estimate the in-hand pose,
correct the erroneous grasps, and successfully insert the component.

Several methods have been developed for pin inspection [8, 9, 10], however,
these methods focus on post-insertion inspection. Pre-insertion inspection is
difficult because of the unknown object position and unknown background. Our
method uses grasping to bring the object to a position with a known
background, and the pose estimation gives the precise position of the pins.
The system can then inspect each pin individually, and is able to reject
components with bent pins successfully.

Another important aspect of a pose estimation systems is the set-up time. The
set-up of a pose estimation system is often a manual task, which is very time-
consuming [11].

Our system introduces a completely automatic set-up procedure based on the CAD
model and grasp pose. Deep learning is used to segment pins, which are used to
find the best configuration to inspect the pins. Based on this configuration,
templates for pose estimation are generated. This procedure gives a set-up
time for a new object as less than one minute.

This paper presents the following main contributions.

  * •

In-hand pose estimation in 2D using stable grasp pose.

  * •

An inspection method to verify pins are straight.

  * •

Deep learning segmentation of object pins.

  * •

Automatic set-up based on CAD model and grasp pose.

  * •

A system for in-hand pose estimation and pin inspection.

The remaining paper is structured as follows: We first review related papers
in Sec. II. In Sec. III, all parts of the developed method are explained. In
Sec. IV, each part is tested individually, and a final insertion system test
is performed. Finally, in Sec. V, a conclusion is given to this paper, and
further work is discussed.

##  II Related Work

This paper presents a novel system for in-hand pose estimation and pin
inspection. However, existing methods also approach the same challenges posed
by this application.

A realm of different methods have been developed to obtain precise pose
estimation, and many methods have also been presented for pin inspection. In
the following, these methods will be detailed, and the reasons that these
methods do not solve the specific challenges addressed by the proposed system
are discussed.

###  II-A Pose Estimation

One method to obtain accurate pose estimation is mechanical fixtures. Here
very precise fixtures are created, and the objects are placed manually to
ensure correct positions. The problem with this solution is that fixtures must
be created for each object used, meaning the fixtures must either be refilled
often or take up large amounts of space. Additionally, manual labour is needed
to fill the fixtures.

A more autonomous solution is mechanical feeders [1]. Here, a large number of
objects are inserted and returned in the correct pose by mechanical design. A
different approach to ensure that the object pose is known is to create
fingers that compensate for errors [2]. The geometry of these fingers guides
the object into a fixed position, thus removing minor errors. This approach
can be combined with the vision to optimize the fingers to specifically
optimize for the vision errors [3]. The disadvantage of these mechanical
solutions is that new hardware must be created for new objects. This severely
limits the flexibility of the systems.

Visual solutions also exist to obtain precise poses. In visual servoing [12],
the robot is continuously updating the position so that the grasped object is
in the desired pose. This is usually performed during the insertion. However,
this does not allow us to perform the pin inspection. Such methods also rely
on training data for the set-up, increasing the set-up time. The object pose
can also be obtained by tracking the object before and after grasping [13].
However, for this approach to work, the object must continuously be observed,
which is hindered by the in-hand occlusion, and the obtained precision is not
high enough for our application.

Another method to predict the in-hand pose is using physical contact [14].
Here the hand is moved towards the table until the object collides. The
position of the collision is used to determine the pose using a particle
filter. This approach has also been extended with a visual check [15].
However, this approach is not feasible for our system as the collisions with
the table can bend the pins.

Poses can also be obtained directly with visual pose estimation. Methods for
6D pose estimation have shown very good results in recent works [16, 5].
However, these methods have challenges with in-hand pose estimation as the
occlusion by the fingers are not accounted for in the method. Additionally,
these methods often lack the very high precision needed for the insertion in
our use case [6].

To solve this, methods have been developed specifically for in-hand pose
estimation [4] [7]. Similar to our method, the hand is removed from the
detection. However, these methods handle larger objects than those in our use
case. They are based on depth information for the pose estimation, which is
difficult to obtain for very small objects.

An approach similar to ours place the object on a table, this table is then
used as a constraint [6]. Like our method, search is restricted to 3D as
constraints related to the object’s movement are used. However, dissimilar to
our method, the used constraint is the position on the table, where our method
uses the constraint of the grasp. This gives our system greater flexibility as
the object can be placed in front of the camera for pin inspection. The pose
estimation developed for this paper is based on a previous method [17]. Here
templates are generated based on the object’s CAD model and matched in the
image using gradients [18]. However, this method performs 6D pose estimation
without constraints, giving uncertainty in the pose estimation, especially in
the depth dimension [19]. Our method can obtain a higher precision as movement
is allowed in fewer dimensions.

###  II-B Pin Inspection

While in-hand pin inspection is a novel approach, pin inspection is a classic
problem. Inspection is generally performed post-insertion on Printed Circuit
Boards (PCB) [8]. Several methods have been developed with cameras placed
above the board [20]. One approach is to extract edges using the Canny Edge
detector [21], and then measure the position [9]. This measurement is similar
to one of the approaches we test in this paper. However, we also test the
match with edge gradients [18]. A different approach to ours is pin inspection
using neural networks [10]. Here a object specific network is trained to
classify if pins are bent or not. We do not train object specific neural
networks in our approach, as training data is not available. As opposed to
these methods our inspection is performed pre-insertion.

Figure 2: Overview of the developed system. Yellow boxes link to the chapter
discussing the respective component in more detail.

##  III Method

The developed system is able to correct for uncertainties in a grasp and
inspect if pins are straight for insertion. The full pipeline consists of two
different parts. An offline phase for set-up of the system and an online pose
estimation and pin inspection phase. Fig. 2 gives an overview of the steps and
the dataflow in the two phases.

The off-line phase requires a CAD model and an approximate grasp pose relative
to the object’s CAD model. The first step is instance segmentation of the
pins. This is performed with a neural network trained on a novel dataset
created for pin segmentation. Based on the pin segmentation and grasp poses
different inspection poses are checked and the best pose is used to generate
templates for the pose estimation.

In the online phase, a pose estimation is first performed to find the object.
Using the pose estimation each individual pin is inspected. If all pins are
determined to be straight the object can be inserted based on the pose
estimation.

###  III-A Pose Estimation

Figure 3:  Visualization of the stability of the in-hand object grasp. The
stable pose of the grasp limits the transform of the object, to translation in
X and Y, and rotation in Z, shown respectively by red, green and blue arrows.

Pose estimation is performed to correct for offsets introduced while grasping
the object. To perform the pose estimation, an image is captured where the
robot has moved the TCP in front of the camera. The pose estimation used in
this work is based on an existing method for 6D pose estimation [17]. This
method is a template matching approach, based on object CAD model. Edge
gradients are matched in the image, and a pyramid search is performed to
reduce the run-time. Finally, pose optimization is performed to refine the
translation.

However, 6D pose estimation with 2D template matching has shown errors as the
depth is unknown [19]. Using stable poses such as table constraints has been
shown to improve precision drastically for 2D template matching [6].

While a table constraint cannot be used in this application, the finger grasp
itself creates a stable pose. The geometry of the fingers and the objects
align the object in the fingers. A planar grasp surface is created which is
oriented in the Z-direction and as such restricts the translation and rotation
of the object. As shown in Fig. 3, the in-hand object is restricted to
translation in X and Y and rotation in the Z-direction. The templates for
matching are created using these constraints, and this allows for a much
smaller search space. This smaller search space improves the precision as
wrong detections cannot be made in the restricted directions. The run-time is
also decreased as fewer templates are searched.

Additionally, as the position of the fingers grasping the object is known,
this ”image area” can be excluded from the template model to increase the
precision.

###  III-B Pin Segmentation

For the full method to work, each pin should be inspected individually. An
instance segmentation of the pins is therefore necessary. While this could be
performed manually, this is a time-consuming and challenging process for an
untrained user, and errors could be introduced by mislabeling pins. The
segmentation of pins is therefore performed automatically. The process is
split into a segmentation using deep learning, and then an instance
segmentation performed with the geometry.

####  III-B1 Network Structure

The object’s CAD model is converted to a point cloud, as several methods for
processing 3D data have shown promising performance on point clouds [22]. To
include details of the small pins, very dense and thereby large point clouds
of the object are needed. The pin processing is, therefore, performed using a
network that has shown good performance for very large point clouds [23]. The
network has also shown an ability to focus on small details, which fits well
with the small pins in our application. As the network in [23] was created for
classification, the segmentation layer from PointNet [24] is added to enable
instance segmentation. The training parameters are identical to the original
paper [23], without the use of error weighing.

####  III-B2 Training Data

A dataset of components with labelled pins has been created to train the
network. The components are downloaded from the free objects from
pcb-3d111https://www.pcb-3d.com/membership_type/free/, using only objects
labelled ”through-hole”. The pin labelling was performed manually, with each
pin segmented. The dataset consists of 895 objects, with an 85/15
training/test split. With 21 % of the points labelled pin. The pin segments
and code for generating the data are available online at
222https://github.com/fhagelskjaer/pin-segmentation.

####  III-B3 Processing CAD model

When processing a CAD model it is first converted to a point cloud and
segmentations are computed using the neural network. The segmentations in the
point cloud are then fit to the original CAD model. For each vertex in the CAD
model a kNN with $k=15$ is used to find the nearest points in the point cloud.
Fifteen points are chosen to smooth the prediction and compensate for errors.
On the CAD model, connected components are then used to compute instance
segmentations of each pin. To ensure that only pins used for insertion are
inspected, pins should be in the positive Z-direction of the TCP. The
resulting segmentation procedure for object D-Sub is shown in Fig. 4.

(a) CAD model.

(b) Point Cloud Representation.

(c) Pin Segmentation.

(d) Instance Segmentation.

Figure 4:  Point cloud segmentation and final instance segmentation of pins.

###  III-C Generating Pose Checks

The object is moved in front of the camera to perform the pose estimation and
pin inspection. The object’s position in front of the camera is known, as the
robot is calibrated to the camera, and the grasp pose is known.

For the pin inspection to be successful, all individual pins should be visible
for the camera. To select the best inspection pose, an occlusion test of all
pins is performed. This test is performed in simulation with the instance
segmentation obtained in Sec. III-B. The object is visualized, and the
visibility is measured for each pin. The visibility is measured by first
projecting the individual pin into the scene and then projecting the entire
object and measuring the occlusion of the pin. To be able to use background
segmentation for the pin inspection, the pin should also not occlude the
object.

As the object can move in-hand, the visibility is checked while rotating the
object from -20 to 20 degrees. The minimum visibility when rotating the object
is then used. Six pre-configured inspection poses are tested this way, and the
pose with the least occlusion is used. Four different inspection poses are
visualized for the D-Sub object in Fig. 5.

(a) Occluded pins.

(b) Pins occluding other pins.

(c) Pins occluding object.

(d) No pin occlusion.

Figure 5:  Pin occlusion test of pre-configured inspection poses.

###  III-D Pin Inspection

To ensure that the object can be inserted, each individual pin needs to be
straight. The pin inspection is based on the pose estimation and the instance
segmentation. Using the pose estimation each pin is inspected in the image,
and the score of the most bent pin is returned. If this score is too low the
component cannot be inserted, and a new component is grasped.

Four different pin inspection methods have been developed and are tested on a
dataset to find the best method. The four methods are gradient check,
intensity overlay, edge overlay, and distance to edge. All methods return a
score and using the test dataset, a cutoff value is found.

For the gradient check, the same check is performed as in the pose estimation
[17]. In the intensity overlay test, a background segmentation is performed,
and the amount of pin overlaying the background is measured. For the edge
overlay edges are computed in the image and the pin edges overlap is measured.
For the distance to edge the image edges are computed and the distance to pin
edges are found. The percentage of edges with less than 2 pixels distance is
then returned.

###  III-E Insertion of objects

The final goal of the system is the insertion of components. To insert the
object the _inserted object_ position is first found. By inserting the object
into the board, and grasping it with the robot the _inserted object_ pose $\
{}_{base}T^{tcp-ins}$ is found. By then moving the robot in front of the
camera and obtaining $\ {}_{cam}T^{obj}$, the pose $\ {}_{tcp}T^{obj}$ is
found as in Eq. 1.

$\ _{tcp}T^{obj}=(\ _{cam}T^{tcp})^{-1}\ _{cam}T^{obj}$ (1)

From this, the position of the _inserted object_ in the base frame can be
computed.

$\ _{base}T^{obj-ins}=\ _{base}T^{ins}\ _{tcp}T^{obj}$ (2)

When an object is grasped and an insertion should be performed the object is
moved in front of the camera and pose estimation is performed. Using the pose
estimation the pins are inspected and objects with bent pins are rejected. If
the object is accepted, the in-hand position, $\ {}_{tcp}T^{obj-current}$, is
found according to Eq. 1. And the new insertion pose can then be calculated as
in Eq. 3.

$\ _{base}T^{tcp-ins-comp}=\\\ \qquad\qquad\ {}_{base}T^{obj-ins}(\
_{tcp}T^{obj-current})^{-1}\\\ $ (3)

##  IV Experiments

To verify the validity of the system, several tests are performed. For each
system component, a test has been performed to measure the individual
validity. Furthermore, a final system test which grasps the object, perform
pose estimation and pin inspection and then inserts the object is reported on.

###  IV-A Pin Segmentation

To process the pin both position and normal-vector information is used for
each point in the point cloud. The point clouds are sampled at a size of 8192
points, and the position information is normalized. The pin segmentation
network is trained for 200 epochs, with a learning rate of 0.001. Two per cent
Gaussian noise is attached to both position and normal vector.

The final model has an accuracy of 96.5 % on the training data with
augmentation, and a test accuracy without augmentation at 97.3 %. The
individual precision and recall for both object and pin is shown in Tab. I. It
is seen that the model obtains a precision of 90 % and 94.8 % recall. In the
training dataset pins represent 17.8 % of the points. To compensate for this,
loss weighing [25] could be used, but the results are adequate for the use
with k-NN.

By using k-NN the small amount of point classification errors are removed when
segmenting the CAD model. The test objects are all segmented correctly.
Segments for the three test objects are shown in Fig. 6.

TABLE I:  Precision-Recall for pin segmentation on the test data.  | Precision | Recall
---|---|---
Component | 98.9 | 97.8
Pin | 90.0 | 94.8

Figure 6:  Point cloud segmentation and final instance segmentation of pins
for the three objects.

###  IV-B Pose Generation

The pose generation step has been introduced to avoid manual set-up, and
thereby have a set-up time that enables fast change-overs. When a new CAD
model is introduced the system should quickly adapt to the object. All tests,
reported on in this subsection, are performed using a Laptop Environment with
an Intel Core i7-10610U CPU, with no GPU.

First, the segmentation is performed using the neural network. The pose checks
are then performed and finally, the templates for the pose estimation are
generated. The run-time of each part for the three objects is shown in Tab.
II.

For our test objects, the set-up time for a new object is less than one
minute. The set-up is, therefore, feasible for fast changeovers in industry.
The largest amount of time is spent on the six inspection poses, thus,
depending on the set-up time requirements, the number of inspection poses
should be limited.

TABLE II:  Set-up time for each object.  Object | Segmentation | Pose Check | Template gen.
---|---|---|---
LED | 6,059s | 11,572s | 2.949s
Pinheader | 6,005s | 29,856s | 7.342s
D-Sub | 6,430s | 40,280s | 7.608s

###  IV-C Pose Estimation

To test the precision of the pose estimation a dataset has been created. The
dataset consist of each object moved to twelve known offset positions along
with a calibration pose. The robot is moved to the calibration pose, $\
{}_{cam}T^{tcp}$, and a pose estimate, ${}_{cam}T^{obj}$, is found. With these
Eq. 4 is used to find the in-hand position of the object.

$\ _{tcp}T^{obj-calib}=(\ _{cam}T^{tcp})^{-1}\ _{cam}T^{obj}$ (4)

At each new position, the $\ {}_{tcp}T^{obj-calib}$ transform is used to find
the expected position in the camera as according to Eq. 5.

$\ _{cam}T^{obj-expected}=\ _{cam}T^{tcp-test}\ _{tcp}T^{obj-calib}$ (5)

By comparing $\ {}_{cam}T^{obj-expected}$ with a new pose estimate at the
position, the offset in position can be found. An illustration of the
different poses and the resulting pose estimation is shown in Fig. 7. The
resulting offset between the expected pose and the pose estimations is shown
in Tab. III. For all objects, the errors is at all times less than one
millimetre.

(a) Calibration pose.

(b) -10 degree angle and 2.5mm offset in X and Y.

(c) 0 degree angle and 2.5mm offset in X and Y.

(d) 10 degree angle and 2.5mm offset in X and Y.

Figure 7: Pose estimation performed on the test dataset. TABLE III:  Error for pose estimation  Object | Mean (mm) | Max (mm)
---|---|---
LED | 0.326 | 0.575
Pinheader | 0.260 | 0.456
D-Sub | 0.267 | 0.529

###  IV-D Pin Inspection

To test the pin inspection a dataset of components with bent and straight pins
has been generated. Different pins have been bent to a degree that the object
cannot be inserted. An example of images from the dataset is shown in Fig. 8.

(a) Straight pins.

(b) Bent pin.

Figure 8:  Visualization of straight and bent pins for the Pinheader.

The resulting score for each of the different pin inspection methods is shown
in Fig. 9. From the scores, the two methods intensity overlay and edge overlay
do not separate bent and straight pins well. The distance to edge works with
Pinheader and D-Sub, but fails with LED. The gradient check method is able to
separate bent and straight pins, but the cutoff value is not similar for the
different objects.

This could stem from inaccuracies in the CAD models. A custom cutoff value is
currently necessary for each object.

Figure 9:  The different pin inspection scores for the three objects. Bent
pins are shown as red crosses and straight pins are shown as green circles.

###  IV-E Complete system test

To test the full usability of the system an insertion test is performed. For
each of the three components, ten objects are tested. Five of the objects have
bent pins and should be rejected by the system. The remaining five should be
inserted successfully.

For all objects with straight pins, the insertion was successful. The bent
components were all rejected by the pin inspection. With the maximum value for
a bent pin and minimum value for a straight pin at LED: 0.82/0.96, Pinheader:
0.49/0.82 and D-Sub: 0.0/0.87.

To further illustrate the system a visual correction is also performed. By
computing ${}_{tcp}T^{obj-current}$ as in Eq. 1, the robot pose to the object
can be placed in the position from the recorded insertion pose, as in Eq. 6.

$\ _{base}T^{tcp-comp}=\ _{base}T^{obj}(\ _{tcp}T^{obj-current})^{-1}$ (6)

Furthermore, as the pin-hole sizes are known the pin-holes can be projected
into the image using Eq. 7. It can then be visually verified that the pose
correction and pin inspection will allow for the insertion of the object.

$x=X\frac{f}{Z}$ (7)

The objects are placed 0.15 meters from the camera with a focal length of 1733
pixels. As the respective hole sizes are 0.9mm, 1.2mm, 1.1mm, for the LED,
Pinheader, and D-Sub, the resulting pixel sizes are 10, 14, and 13. In Fig. 10
and Fig. 11 the visual corrections are shown the Pinheader and D-Sub,
respectively.

(a) Initial Grasp.

(b) Corrected Pose.

(c) Initial Grasp.

(d) Corrected Pose.

(e) Initial Grasp.

(f) Corrected Pose.

Figure 10:  Visualization of the pose correction and insertion for object
Pinheader. The holes are projected into the scene based on the insertion pose
and coloured green. Holes colliding with the pins are coloured red. Best
viewed digitally.

(a) Initial Grasp.

(b) Corrected Pose.

(c) Initial Grasp.

(d) Corrected Pose.

Figure 11:  Visualization of the pose correction and insertion for object
D-Sub. The holes are projected into the scene based on the insertion pose and
coloured green. Holes colliding with the pins are coloured red. Best viewed
digitally.

##  V Conclusion

In this work, we presented a novel system for in-hand pose estimation and pin
inspection. The system consists of several novel contributions. A network
trained for segmentation of object pins. A method to evaluate inspection poses
based on robot configurations. In-hand pose estimation based on stable grasp
poses. And a viable pin inspection method. Finally, the parts are connected to
perform a successful insertion task.

For further work, the system could be extended with multiple cameras to
inspect pins from multiple views. The pin inspection could also be performed
by projecting the expected holes into the scene and checking if the object can
be inserted. Furthermore, actual pose estimation of the pins could be
performed, and the insertion strategy could be based on the pin positions as
opposed to the object position. This would also allow for the development of
more advanced insertion strategies that could potentially correct for
significantly bent pins.

The in-hand pose estimation could be tested on other objects, and extended to
work with multiple stable grasp poses.

## Acknowledgement

This project was funded in part by Innovation Fund Denmark through the project
MADE FAST, in part by the SDU I4.0-Lab.

## References

  * [1] S. Mathiesen, L. C. Sørensen, D. Kraft, and L.-P. Ellekilde, “Optimisation of trap design for vibratory bowl feeders,” in _2018 IEEE international conference on robotics and automation (ICRA)_. IEEE, 2018, pp. 3467–3474.
  * [2] A. Wolniakowski, A. Gams, A. Kramberger, and A. Ude, “Compensating pose uncertainties through appropriate gripper finger cutouts,” _acta mechanica et automatica_ , vol. 1, no. 12, pp. 78–83, 2018.
  * [3] F. Hagelskjær, A. Kramberger, A. Wolniakowski, T. R. Savarimuthu, and N. Krüger, “Combined optimization of gripper finger design and pose estimation processes for advanced industrial assembly,” in _2019 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2019, pp. 2022–2029.
  * [4] C. Choi, J. D. Preto, and D. Rus, “Using vision for pre-and post-grasping object localization for soft hands,” in _International Symposium on Experimental Robotics_. Springer, 2016, pp. 601–612.
  * [5] T. Hodan, M. Sundermeyer, B. Drost, Y. Labbe, E. Brachmann, F. Michel, C. Rother, and J. Matas, “Bop challenge 2020 on 6d object localization,” _arXiv preprint arXiv:2009.07378_ , 2020.
  * [6] F. Hagelskjær, T. R. Savarimuthu, N. Krüger, and A. G. Buch, “Using spatial constraints for fast set-up of precise pose estimation in an industrial setting,” in _2019 IEEE 15th International Conference on Automation Science and Engineering (CASE)_. IEEE, 2019, pp. 1308–1314.
  * [7] B. Wen, C. Mitash, S. Soorian, A. Kimmel, A. Sintov, and K. E. Bekris, “Robust, occlusion-aware pose estimation for objects grasped by adaptive hands,” in _2020 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2020, pp. 6210–6217.
  * [8] A. Kashitani, N. Takanashi, and N. Tagawa, “A solder joint inspection system for surface mounted pin grid arrays,” in _Proceedings of IECON’93-19th Annual Conference of IEEE Industrial Electronics_. IEEE, 1993, pp. 1865–1870.
  * [9] Y. Huang, Q. Pan, Q. Liu, X. P. He, Y. F. Liu, and Y. Q. Yu, “Application of improved canny algorithm on the ic chip pin inspection,” in _Advanced Materials Research_ , vol. 317. Trans Tech Publ, 2011, pp. 854–858.
  * [10] I. Edinbarough, R. Balderas, and S. Bose, “A vision and robot based on-line inspection monitoring system for electronic manufacturing,” _Computers in Industry_ , vol. 56, no. 8-9, pp. 986–996, 2005.
  * [11] F. Hagelskjær, A. G. Buch, and N. Krüger, “Does vision work well enough for industry?” in _Proceedings of the 13th International Joint Conference on Computer Vision, Imaging and Computer Graphics Theory and Applications_ , vol. 4. SCITEPRESS Digital Library, 2019, pp. 198–205.
  * [12] C. Yu, Z. Cai, H. Pham, and Q.-C. Pham, “Siamese convolutional neural network for sub-millimeter-accurate camera pose estimation and visual servoing,” in _2019 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_. IEEE, 2019, pp. 935–941.
  * [13] J. Liang, A. Handa, K. Van Wyk, V. Makoviychuk, O. Kroemer, and D. Fox, “In-hand object pose tracking via contact feedback and gpu-accelerated robotic simulation,” in _2020 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2020, pp. 6203–6209.
  * [14] F. Von Drigalski, S. Taniguchi, R. Lee, T. Matsubara, M. Hamaya, K. Tanaka, and Y. Ijiri, “Contact-based in-hand pose estimation using bayesian state estimation and particle filtering,” in _2020 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2020, pp. 7294–7299.
  * [15] F. von Drigalski, K. Hayashi, Y. Huang, R. Yonetani, M. Hamaya, K. Tanaka, and Y. Ijiri, “Precise multi-modal in-hand pose estimation using low-precision sensors for robotic assembly,” in _2021 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2021, pp. 968–974.
  * [16] F. Hagelskjaer and A. G. Buch, “Parapose: Parameter and domain randomization optimization for pose estimation using synthetic data,” _arXiv preprint arXiv:2203.00945_ , 2022.
  * [17] M. Ulrich, C. Wiedemann, and C. Steger, “Combining scale-space and similarity-based aspect graphs for fast 3d object recognition,” _IEEE transactions on pattern analysis and machine intelligence_ , vol. 34, no. 10, pp. 1902–1914, 2012.
  * [18] C. Steger, “Occlusion, clutter, and illumination invariant object recognition,” _International Archives of Photogrammetry Remote Sensing and Spatial Information Sciences_ , vol. 34, no. 3/A, pp. 345–350, 2002.
  * [19] B. Drost, M. Ulrich, P. Bergmann, P. Härtinger, and C. Steger, “Introducing mvtec itodd-a dataset for 3d object recognition in industry.” in _ICCV Workshops_ , 2017, pp. 2200–2208.
  * [20] Á. Becker and L. Jakab, “Optical examination of pin insertion process,” _Periodica Polytechnica Electrical Engineering (Archives)_ , vol. 53, no. 1-2, pp. 85–89, 2009.
  * [21] J. Canny, “A computational approach to edge detection,” _IEEE Transactions on pattern analysis and machine intelligence_ , no. 6, pp. 679–698, 1986.
  * [22] Y. Guo, H. Wang, Q. Hu, H. Liu, L. Liu, and M. Bennamoun, “Deep learning for 3d point clouds: A survey,” _IEEE transactions on pattern analysis and machine intelligence_ , 2020.
  * [23] F. Hagelskjaer, “Deep learning classification of large-scale point clouds: A case study on cuneiform tablets,” _arXiv preprint arXiv:2202.10851_ , 2022\.
  * [24] C. R. Qi, H. Su, K. Mo, and L. J. Guibas, “Pointnet: Deep learning on point sets for 3d classification and segmentation,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2017, pp. 652–660.
  * [25] Y.-X. Wang, D. Ramanan, and M. Hebert, “Learning to model the tail,” in _Proceedings of the 31st International Conference on Neural Information Processing Systems_ , 2017, pp. 7032–7042.


# Control of energy spectra and enhancement of energy conversion of fast
electrons generated by dual-color picosecond lasers

Tie-Huai Zhang  Beijing National Laboratory for Condensed Matter Physics,
Institute of Physics, CAS, Beijing 100190, China  School of Physical Sciences,
University of Chinese Academy of Sciences, Beijing 100049, China  Wei-Min Wang
<EMAIL_ADDRESS>Department of Physics and Beijing Key Laboratory of
Opto-electronic Functional Materials and Micro-nano Devices, Renmin University
of China, Beijing 100872, China  Key Laboratory of Quantum State Construction
and Manipulation (Ministry of Education), Renmin University of China, Beijing,
100872, China  IFSA Collaborative Innovation Center, Shanghai Jiao Tong
University, Shanghai 200240, China  Yu-Tong Li<EMAIL_ADDRESS>Beijing
National Laboratory for Condensed Matter Physics, Institute of Physics, CAS,
Beijing 100190, China  School of Physical Sciences, University of Chinese
Academy of Sciences, Beijing 100049, China  IFSA Collaborative Innovation
Center, Shanghai Jiao Tong University, Shanghai 200240, China  Songshan Lake
Materials Laboratory, Dongguan, Guangdong 523808, China  Jie Zhang  Beijing
National Laboratory for Condensed Matter Physics, Institute of Physics, CAS,
Beijing 100190, China  IFSA Collaborative Innovation Center, Shanghai Jiao
Tong University, Shanghai 200240, China  Key Laboratory for Laser Plasmas
(MoE) and School of Physics and Astronomy, Shanghai Jiao Tong University,
Shanghai 200240, China

###### Abstract

After the successful fusion ignition at National Ignition Facility, seeking
for a high-gain fusion scheme becomes the next hot-spot in inertial
confinement fusion community. Fast ignition provides an alternative due to its
potential to reduce the energy of driven lasers and achieve higher target
gain, whose core step is generating fast electron beam using picosecond
lasers. The properties of the electron beam, such as energy spectra,
determines the succeed heating process and the yield. In this study, with 2D
particle-in-cell simulations, we find that the energy transfer ratio would
significantly increase, and the energy spectra would be lifted up in the dual-
color injecting scheme, where an extra weak low-frequency laser is injected
with the main pulse. These phenomena are attributed to the longitudinal
electric field modulation in low-density region and electromagnetically-
induced-transparency process in the microchannels formed near the critical
surface. Our results can be applied to fast ignition schemes and experimental
designs of super-hot electron generation with different pulse durations. For
example, in the DCI project with a fast ignition scheme, such a dual-color
intense picosecond laser facility is under construction to control the fast
electron generation and enhance the laser transparency.

##  I Introduction

Since 2021, several breakthroughs have been achieved in inertial confinement
fusion (ICF) experiments at National Ignition Facility, marking that the
threshold of controlled fusion ignition is approached for the first time
Zylstra  _et al._ (2022); Abu-Shawareb  _et al._ (2022). However, the extra
energy transfer process at the hohlraum would lead to substantial radiation
loss, limiting its steps towards high target gain. Furthermore, there are
still problems such as energy conversion efficiency and the balance on
compression and hydrodynamic instability Colaïtis  _et al._ (2022);
Schlossberg  _et al._ (2021). Besides central ignition, fast ignition (FI)
Tabak  _et al._ (1994); Kodama  _et al._ (2001); Tabak  _et al._ (2006);
Kitagawa  _et al._ (2022) is another possible scheme for ICF, which drives the
deuterium–tritium (DT) fuel with weaker nanosecond laser beams compared to the
central ignition schemes and uses picosecond pulse(s) to generate super-hot
electron beam and ignite the compressed plasma. In this scheme, the energy
conversion from laser to super-hot electrons lies in the central position of
the whole experiment design. To shorten the transport distance of electron
beam and reduce energy losses, a gold cone is widely used in many FI schemes
such as cone-in-shell Nagatomo  _et al._ (2019) and double-cone ignition (DCI)
Zhang  _et al._ (2020) to keep the ps pulse apart from low-density corona
plasma, guiding laser pulse towards the core Nagatomo  _et al._ (2019); Rusby
_et al._ (2021).

The energy spectrum of super-hot electron beam is another key factor in FI
scheme, since the electrons have to travel a long distance before reaching the
core plasma. The position of energy deposition and hot spot generation are
very sensitive to the electron spectrum, which needs detailed optimization for
specific schemes. In typical cone-in-shell scheme, this transport distance is
around 80 µm Robinson  _et al._ (2014); Nagatomo  _et al._ (2021), which
implies that only a part of super-hot electrons lying in the specific energy
range can heat the hot-spot region effectively.

Apart from ICF, the electron beam generated by intense short-pulse (sub-
picosecond or femtosecond) laser interacting with solid target also has a wide
applications such as generating broadband x-ray emission Sawada  _et al._
(2020) and strong terahertz radiation Woldegeorgis  _et al._ (2019). There are
various temperature scaling laws of super-hot electron generated by laser-
plasma interaction (LPI), distinguished from pulse duration and laser
intensity Gibbon and Bell (1992); Forslund  _et al._ (1977); Wilks  _et al._
(1992). When laser intensity reaches the relativistic region, $\bf J\times B$
heating is the dominant acceleration mechanism and the ”ponderomotive
temperature” [Eq. (1)] is often used to estimate electron temperature
generated by a linear-polarized laser Wilks  _et al._ (1992), where $I_{18}$
is the intensity in $10^{18}$ W/cm2.

$T_{p}\simeq 0.511\left[\sqrt{1+0.73(I_{18}\lambda_{\mu
m}^{2})}-1\right]\quad(\text{MeV})$ (1)

However, the pulse duration and density scale length would strongly affect
temperature of super-hot electrons Ivanov  _et al._ (2014); Yabuuchi  _et al._
(2010). Actually, the electron spectrum depends on much more parameters
besides laser intensity and wavelength, and the electron energy does not obey
the single-temperature Maxwell distribution. To improve the yields of
electrons in aiming energy range, further investigations are urgently
demanded. Hence, various strategies have been raised to modify the structure
of super-hot electron spectra and improve the energy coupling efficiency.

One of the strategies to improve optical injection of super-hot electrons is
introducing the second laser pulse. For example, a counter-propagating
additional laser pulse can result in an increase on injected electron energy
and number Wang and Sheng (2008), or excite beat wave for generating mono-
energetic electron beam Faure  _et al._ (2006). An orthogonal beam can also
dephase the electrons trapped in the wake-field and inspire transverse
emittance in the accelerated bunch Dodd  _et al._ (2004). For solid target,
limited by the laser arrangement, the coaxial schemes Chitgar  _et al._
(2020); Zeng  _et al._ (2015); Yu  _et al._ (2014) and small-angle
configurations Morace  _et al._ (2019); Zhou  _et al._ (2022) are more
applicable. To improve the electron beam, the short-wavelength pulse, whose
dimensionless field strength is smaller with the same intensity ($a\propto
I\lambda^{-2}$), is also widely used to generate a softer energy spectrum. The
short-wavelength pulse, corresponding to a higher critical density, can also
be applied to cavitate ponderomotive channels in over-dense region when the
pulse duration is long enough, guiding the long-wavelength pulse into these
micro structures and obtaining a more sufficient LPI. Moreover, due to the
lower diffraction limit, the waist of laser can be narrowed further, reaching
a higher concentration of laser energy Scott  _et al._ (2013). Apart from
that, using short-wavelength pulse to generate super-hot electrons in FI can
reduce the possible negative effect from the preplasma in gold cone Baton  _et
al._ (2008).

By combining the two strategies above, i.e., using dual-color lasers, it is
possible to further enhance the quality of electron beams. For example, a
coaxial long-wavelength pulse followed by a short-wavelength pulse can produce
a quasi-monoenergetic electron beam with low emittance Yu  _et al._ (2014),
and a femtosecond laser with fundamental frequency and its harmonic co-
transmit in mixed gas plasma can trigger period electron injection, generating
multichromatic narrow energy-spread electron brunches Zeng  _et al._ (2015).

A common method to generate a short-wavelength pulse is frequency-doubling.
The conversion efficiency of nonlinear crystals can reaches $\sim 65\%$, and
the rest energy remains as fundamental frequency component. Therefore, whether
we should retain or filter out the fundamental frequency component after
frequency-doubling in fast igintion, such as the DCI project in China, is an
urgent realistic question.

In this paper, we investigated the effects of simultaneous dual-color laser
injection on solid thick targets and the underlying physical processes. In
Section II, the simulation parameters and corresponding electron spectrum and
energy transfer efficiency are given in detail, showing a significant yield
increase of high-energy range electrons in dual-color injecting cases. Section
III and IV analyze the possible mechanisms of dual-color effects with
different pulse duration. Finally, the polarization and frequency relationship
are discussed in Section V, followed by the conclusion.

##  II 2D Particle-in-Cell Simulation

We assume that the maximum approachable intensity is $I_{0}=5\times 10^{19}$
W/cm2, and the 2nd harmonic ($2\omega$) conversion efficiency is 65% (for
typical nonlinear crystal like BBO, this value can reach up to 60—70%). Part
of fundamental frequency component is retained here to investigate the
influence on laser-plasma coupling and electron energy spectra in dual-color
injecting cases.

###  II.1 Simulation Parameters

2D PIC simulations are performed with the KLAPS codeWang  _et al._ (2015). The
laser beams in our simulations transmit along $+x$ axis. Every beam has an
intensity profile of $I(y,t)=I_{\alpha}e^{-y^{2}/\sigma^{2}}\sin^{2}(\pi
t/\tau),\alpha=1,2$, with p-polarization unless otherwise specified. The
footnote $\alpha$ denotes parameters of different laser beams in a single
simulation case. In this work, we specify Laser 1 as the fundamental frequency
beam, where $\lambda_{1}=1$ µm, and Laser 2 as the $2\omega$-beam. All the
related parameters, for instance, the critical density
$n_{cr}=\varepsilon_{0}m_{e}\omega^{2}/c^{2}$ and dimensionless strength ${\bf
a}=e{\bf E}/m_{e}c\omega$, are calculated with the frequency of the first
pulse. The target is set to be collisionless Au plasma with a static
ionization of $Z=10$. The initial density profile $n_{0}(x,y)=n_{0}(x)$ is
shown as Eq. (2), which has an exponential pre-plasma due to the ASE pedestal.

$\displaystyle n_{0}(x)=\left\\{\begin{aligned}
&n_{min}\exp(\frac{x-x_{0}}{L}),&&x_{0}<x<x_{1}\\\ &n_{max},&&x_{1}<x<x_{2}\\\
&0,&&\text{Others}\end{aligned}\right.$ (2)

The maximum density, left plasma-vacuum boundary and laser waist are fixed at
$n_{max}=100n_{c}$, $x_{0}=3$ µm and $\sigma=10$ µm. To hold the continuity of
density, we take $x_{1}=x_{0}+L\ln(n_{max}/n_{min})$.

To systematically investigate the effects of dual-color injection with
different pulse duration, we set $\tau=1\text{ ps }(300T_{0}),0.2\text{ ps
}(60T_{0})$ and $66.7\text{ fs }(20T_{0})$ for different cases, where $T_{0}$
is the field period of the fundamental frequency laser. The $\tau$s above
correspond to typical picosecond, sub-picosecond and femtosecond laser
duration separately. Considering that the prepulse from a long pulse is
usually stronger than a short one, the simulation box of the picosecond case
is $64\times 48$ µm ($2048\times 1536$ cells), with $L=8$ µm,
$n_{min}=0.1n_{c}$, and $x_{2}=63$ µm. For the sub-picosecond and femtosecond
cases, the simulation box is $40\times 32$ µm ($1280\times 1024$ cells), with
$L=4$ µm, $n_{min}=0.01n_{c}$, and $x_{2}=39$ µm. It should be note that scale
length variation cause by the change of laser wavelength Scott  _et al._
(2013) is not considered here.

Detecting surfaces perpendicular to $x$ axis are placed at regular intervals
of $\Delta x=2$ µm to record the energy flux of the super-hot electrons which
pass through the surface with $p_{x}>0,E_{k}>50$ keV for the first time. Then,
the spectra of electrons recorded by the detecting surface near $x_{1}$
(located at $x=57$ µm or $x=31$ µm) are given in Section II.2.

###  II.2 Simulation Results

The super-hot electron energy flux and spectrum in different injecting cases
are shown in Fig. 1, where $I_{1}=I_{\omega},I_{2}=I_{2\omega}$. Figure 1(a)
infers that the dual-color injecting strategy can lift the energy coupling
efficiency. A 6% increase can be found at the detecting surface near
$n_{0}=100n_{c}$ with $5\%I_{0}$ fundamental frequency laser doping. When the
fundamental frequency component is mixed in the $2\omega$-beam, the high-
energy part of electron energy spectrum has a significant rise, while the low-
energy electron number falls down a little. In the case of
$I_{\omega}=30\%I_{0},I_{2\omega}=65\%I_{0}$, the spectrum curve has become
very close to the fundamental frequency case in the high-energy range.
Actually, the two curves intersect at 4.6 MeV and 12.2 MeV, which means that
in both low- and high-energy regions, the dual-color injecting scheme performs
better than fundamental frequency case even though there is an energy loss of
5% in consideration of the frequency-doubling process.

Figure 1:  (a) Super-hot electron flux passing through detecting surfaces in
different injecting cases. The data have been normalized with the incident
total laser energy. (b) The energy spectra of super-hot electrons recorded by
the surface at $x=57$ µm ($n_{0}=85.4n_{c}$), corresponding to the dots marked
in Figure 1(a).

Furthermore, we investigate the electron spectra in cases of shorter pulse
duration, and obtain the similar results (Fig. 2). The curves of dual-color
injecting cases are also lifted in the high-energy range comparing to the pure
$2\omega$ case for both femtosecond and sub-picosecond-laser injection. It
should be pointed out that $I_{\omega}=30\%I_{0},I_{2\omega}=65\%I_{0}$ case
behaves better than or almost the same as fundamental frequency case over the
whole energy range for sub-picosecond-laser injection [Fig. 2(b)].

Figure 2:  The energy spectra of super-hot electrons recorded by the surface
at $x=39$ µm ($n_{0}=81.0n_{c}$), with the laser duration (a) $\tau=66.7$ fs,
(b) $\tau=0.2$ ps

The mechanism of dual-color injecting effect is going to be analyzed in detail
in the following two sections. Nevertheless, a caveat must be addressed in
advance. Although the lower-frequency beam is usually expected to generate
super-hot electrons with higher kinetic energy, the intensity also plays a
part. According to the scaling estimation Eq. (1), the ponderomotive
temperature is related to $I\lambda^{2}$, which means that it is improper to
consider the two beams’ contribution separately and add the electrons produced
by every single laser together linearly. Taking the
$I_{1}=5\%I_{0},I_{2}=65\%I_{0}$ case as a counter-example,
$(I_{1}\lambda^{2}_{1})/(I_{2}\lambda^{2}_{2})=0.31<1$, the fundamental
frequency laser corresponds to a even lower ”ponderomotive temperature”, while
the high-energy part of the spectrum curve soars compare to the $2\omega$-beam
injection case. In a word, investigations of the associative effects of co-
transmission of dual-color beams in plasma are inevitably required.

##  III Longitudinal Electric Field Excited by the Dual Wavelength Laser Field

Since the increase of high-energy part of the electron flux appears in cases
with various pulse duration scale, it is deduced that mixture of the two beams
would alter the electron acceleration process from the linear stage of laser-
plasma interaction. The phase space diagram and dimensionless field strength
$a_{x}$ on $x$ axis of two different femtosecond-laser injecting schemes,
$I_{1}=5\%,I_{2}=65\%I_{0}$ (dual-wavelength beams) and
$I_{1}=0,I_{2}=65\%I_{0}$ (pure $2\omega$ beam), are compared in Fig. 3.

Figure 3:  The phase space diagram ($p_{x}-x$) and dimensionless field
strength $a_{x}$ on $x$ axis of (a), (c) $I_{1}=5\%,I_{2}=65\%I_{0}$ case and
(b), (d) $I_{1}=0,I_{2}=65\%I_{0}$ case at $t=33.3$ fs ($10T_{0}$). The data
have been normalized with the incident total laser energy. Figure (a) and (b)
plot 0.5% of the electrons randomly whose kinetic energy are larger than 50
keV, with a shared coordinate axis. The size and color of points denote the
weight of picked simulation particles.

With a $5\%I_{0}$ mixture of $2\omega$ laser, the phase space diagram deviates
from a modulated structure symmetric about $p_{x}=0$ [Fig. 3(b)]. Several
peaks and sub-peaks emerge with certain space period ($\sim\lambda_{1}$),
which correspond to electrons with much higher forward momentum and kinetic
energy. Changes are also found in longitudinal electric field diagrams.
Sawtooth structure appears on $a_{x}$ plot shown in Fig. 3(c), with a matched
spacial period as peaks in phase space diagram. When the intense laser field
interacts with underdense plasma, a longitudinal wakefield would be exited
behind the laser front, which can trap electrons and accelerate them as the
wave propagates Tajima and Dawson (1979); Brantov  _et al._ (2008). The quasi-
static approximation (QSA) can be applied to plane wave propagating in uniform
underdense cold plasma, until the wave envelope changes significantly, whose
characteristic timescale is usually the Rayleigh diffraction time Gibbon
(2005). In this scheme, there is a coordinate transformation
$\xi=x-v_{p}t,\tau=t$ in 1D problem and the longitudinal field
$a_{x}(\xi)=-\differential{\phi}/\differential{(k_{p}\xi)}$ satisfies
Poisson’s equation Dalla and Lontano (1993):

$\derivative[2]{\phi}{\xi}=k_{p}^{2}\gamma^{2}\left[\frac{\beta_{p}(1+\phi)}{\sqrt{(1+\phi)^{2}-\gamma^{-2}(1+a^{2})}}-1\right],$
(3)

where $\beta_{p}\sim 1$ is the phase velocity of $a_{x}$ excited by the laser
field, $\gamma=1/\sqrt{1-\beta_{p}^{2}}$ and $k_{p}=\omega_{p}/c$.
Particularly, in rarefied plasma ($n\ll n_{c}$), $\beta_{p}=1$ and Eq. (3) can
be reduced to Eq. (4) V.  _et al._ (1989):

$\derivative[2]{\phi}{\xi}=k_{p}^{2}\frac{(1+a^{2})-(1+\phi)^{2}}{2(1+\phi)^{2}}.$
(4)

Once the pulse envelope is substituted into the Poisson’s equation, the well-
known wakefield solution can be obtained numerically V.  _et al._ (1989);
Dalla and Lontano (1993), which is theoretical basis of wakefield
acceleration. Here, since the envelopes of laser beams are the same, the
oscillating term of field must be taken into consideration. For the validity
of QSA, it is assumed that both of the two-color beams have the same phase
velocity $c$, suggesting that the dispersion is ignorable and
$\omega_{1}/\omega_{2}=k_{1}/k_{2}$.

When we normalize $k_{1}\xi\rightarrow\xi$, and take

${\bf a}(\xi)=\left\\{\begin{aligned}
&\sin\left(\frac{\xi}{2N}\right)\left[a_{1}\hat{\bf
e}_{1}\sin(\xi)+a_{2}\hat{\bf
e}_{2}\sin\left(\frac{k_{2}}{k_{1}}\xi\right)\right],&&0<\xi<2\pi N\\\
&0,&&Others\end{aligned}\right.,$ (5)

where $N=\tau/T_{1}$ is the pulse duration normalized by field circulation of
Laser 1, the longitudinal field can be solved numerically. Figure 4 compares
$a_{x}$s of the dual-color (a) and pure-$2\omega$ injecting case given by Eq.
(4). The consideration of full laser field brings fluctuations on the smooth
$a_{x}$ curve (shown as the dashed line in Fig. 4), which matches the PIC
simulation result of Fig. 3(d). In the dual-color case, The sawtooth structure
is observed clearly, leading to the deduction that co-transmission of dual-
color transverse field can modulate the longitudinal static electric field,
forming periodic sawtooth structure. The exited positive longitudinal field
would hamper forward acceleration. However, for electrons located at the
troughs of modulated wave, they suffer a decreased damping force and can gain
a higher forward momentum, corresponding to the peaks in Fig. 3(a). The
modulating strength is stronger with a mixture of higher amplitude fundamental
frequency beam. Therefore, the high-energy part of electron energy spectrum
curve is lifted, resulting in harder spectrum and a higher conversion
efficiency.

In longer duration cases, the mechanism discussed above is also valid.
However, since the intensity grows much slower than ultra-short pulse, these
electrons with relatively higher initial forward momentum at the wave troughs
would be further accelerated by subsequent laser field directly or other
mechanisms like stochastic acceleration Sheng  _et al._ (2002); Meyer-ter Vehn
and Sheng (1999).

Figure 4:  The numerical solutions of Eq. 4 of (a) dual-color case
($a_{1}=2a_{2}=5$), and (b) pure-$2\omega$ case ($a_{2}=0$). Both of the two
cases take $\hat{\bf e}_{1}=\hat{\bf e}_{2}=\hat{\bf e}_{y}$, $k_{2}/k_{1}=2$,
$n_{c}/n=25$, and $N=20$. The solid lines denote that $a$ takes the full
version of Eq. (5), and the dashed lines correspond to solutions derived from
a laser field envelope.

##  IV Long Pulse: Microchannels Forming on the Critical Surface

Although modulation of $a_{x}$ can generate a part of electrons as seed for
later acceleration, the longitudinal field far behind the wavefront, where the
laser amplitude actually rises in long pulse cases, would become chaotic,
leading to the destruction of sawtooth structure. Besides, there is enough
time for plasma to respond the perturbation from laser fields and for
nonlinear effects to develop. One of the direct consequences is that a set of
microchannels is cavitated due to the perturbation of electron density by
ponderomotive force and transverse instabilities Huang  _et al._ (2013); Wan
_et al._ (2018). Laser fields can propagate along these microchannels Sprangle
_et al._ (1992); Pukhov and Meyer-Ter-Vehn (1996), and thus a longer distance
for LPI is approachable, resulting in gains on energy coupling efficiency. The
acceleration of electrons in these channels would also be affected by the
static fields, where some of electrons are expelled out of the channel center,
forming a static electron field $E_{y}^{s}\propto y$ and $B_{z}^{s}\propto-y$
Pukhov  _et al._ (1999); Arefiev  _et al._ (2012).

Figure 5 shows the electron density, normalized electric field $a_{y}$ and
normalized static magnetic field $\langle b_{z}\rangle$ in different
picosecond simulation cases at $t=400$ fs. The magnetic fields are averaged
over two laser field periods. To show the details of laser-plasma interaction
region, the plot spacial range is limited to $X\times
Y\in[15,35]\times[-15,15]$ (µm). It is confirmed with the density diagrams
that the over-dense region where $n_{0}>n_{cr}$ is perturbed, displaying some
certain parallel structures. There are also channel structures in electric
field plots which overlap the perturbed area in density diagrams, implying
that laser fields do propagate along the microchannels. Hence, super-hot
electrons can be generated and transport in these channels towards the rear
surface of the target, which is consistent with the static magnetic field
configuration in Fig. 5(g)-(i). It should be noticed that the width of
microchannels varies among different cases, which could be related to
parameters such as the normalized field strength, plasma density and scale
length. Obviously, a longer and wider channel is beneficial to LPI and energy
coupling, which can be considered as one of the reasons that dual-color
injecting scheme has a higher coupling efficiency.

Figure 5:  Diagrams of (a)-(c) electron density, (d)-(f) normalized electric
field $a_{y}$, and (g)-(i) normalized static magnetic field $\langle
b_{z}\rangle$ in different picosecond simulation cases at $t=400$ fs. The
first row: $I_{1}=100\%I_{0},I_{2}=0$ (pure fundamental frequency); the second
row: $I_{1}=30\%I_{0},I_{2}=65\%I_{0}$ (dual-color); the second row:
$I_{1}=0,I_{2}=65\%I_{0}$ (pure $2\omega$). The black and white dash-dotted
lines in (a)-(c) correspond to $n_{0}(x)=n_{cr}$ and $n_{0}(x)=4n_{cr}$
(critical surface of $2\omega$ beam) separately.

Although plasma in microchannels is perturbed and expelled by the laser field,
the fundamental frequency component can still penetrate into the over-dense
region. In the pure fundamental frequency case, one can ascribe this
phenomenon to the well-known relativistic transparency (RT) effect Kaw and
Dawson (1970). Since the zeroth-order dispersion relationship for relativistic
laser field is
$\omega^{2}=k^{2}c^{2}+\omega_{p}^{2}/\gamma_{0},\gamma_{0}=\sqrt{1+\langle
a^{2}\rangle/2}$, the effective plasma frequency decreases as the rise of
normalized field amplitude, making plasma transparent for the electromagnetic
wave.

While in the dual-color injecting case, even though the intensity of
fundamental frequency component is as low as $5\%I_{0}$, it can be detected
where the local electron density rises up to several times more than
$n_{cr1}$. Figure 6 shows the time-dependent evolution of smoothed $n_{e}$ and
filtered normalized electric field $a_{y}$ at $(x,y)=$(31 µm, 0). Between 500
fs to 800 fs, strong fundamental frequency signals can be detected even if the
local $n_{e}$ is much higher than the RT-modified critical density
$\gamma_{1}n_{cr1}=1.63n_{cr1}$, which can be attributed to the relativistic
electromagnetically induced transparency (EIT) mechanism Tie . In the physical
picture of EIT, when a low-frequency laser $(a_{1},\omega_{1})$ co-transmit
with a pump laser with $\omega_{2}>\omega_{p}$ in over-dense plasma, the
induced current of Laser 1 tends to be canceled by the beat current of pump
wave and plasma wave, leading to the transparency of Laser 1 Harris (1996);
Gordon  _et al._ (2000a, b).

Figure 6:  Smoothed electron density $n_{e}$ and normalized laser fields
$a_{y}$ filtered within $[0.95\omega_{0},1.05\omega_{0}]$ at $(x,y)=$(31 µm,
0) in case $I_{\omega}=5\%I_{0},I_{2\omega}=65\%I_{0}$. Since the cavitation
process would push the plasma forward, the local density rises from
$3.32n_{cr1}$ to $5.78n_{cr1}$ in $750$ fs, then falls down to $2.73n_{cr1}$.

At the late stage of LPI, the microchannels would expand transversely and
break down. Eventually, these channels merge and forming a cavity with the
radius $r\sim\sigma$, which can be converted to the typical hole-boring
picture Wilks  _et al._ (1992); Pukhov and Meyer-Ter-Vehn (1997). However, the
laser intensity is not extremely strong and pulse duration is relatively short
in this paper, the subsequent hole-boring process is weak and can be omitted
here Kemp  _et al._ (2008); Weng  _et al._ (2012).

##  V Dependency on polarization and frequency relationship

In Section II, III and IV, both of the two beams are fixed at p-polarization.
In some cases, the nonlinear crystals can generate second harmonic waves with
the polarization direction perpendicular to the pump beam. To investigate the
effect of beam polarization on electron energy spectrum, we changed the
polarization direction of Laser 1 to $\hat{\bf e}_{1}=\hat{\bf e}_{z}$
(s-polarization) and took the energy proportion
$I_{1}=30\%I_{0},I_{2}=65\%I_{0}$ in dual-color injecting scheme, showing the
electron energy spectrum of femtosecond and picosecond cases in Figure 7. The
orthogonal polarization directions have a negative effect on the super-hot
electron generation in both of the two cases, especially in high-energy
region. It is unsurprising since one can calculate the modulation of
longitudinal field is greatly reduced when $\hat{\bf e}_{1}\cdot\hat{\bf
e}_{2}=0$. According to Eq. (5) and letting $a_{1}=\epsilon a_{2}$,
$a^{2}={\bf a}\cdot{\bf a}=\sin^{2}(k_{2}\xi/k_{1})+(\hat{\bf
e}_{1}\cdot\hat{\bf
e}_{2})\epsilon\sin(\xi)\sin(k_{2}\xi/k_{1})+\epsilon^{2}\sin^{2}(\xi)$, and
$\hat{\bf e}_{1}\cdot\hat{\bf e}_{2}=0$ leads to the first-order perturbation
term being zero. In long-pulse injecting case, the microchannels still exist,
where an external driving term appears in parametric instability equations
Arefiev  _et al._ (2012) due to the orthogonal polarization. Since the induced
currents of the two laser beams fluctuate in different directions, the EIT
effect is expected to vanish Tie , resulting in the decrease of super-hot
electron output. However, it should be noted that the simulations here are in
2D geometry, where the inhomogeneous effects along $z$ direction cannot appear
and need further investigation of 3D simulations.

Figure 7:  The energy spectrum of (a) femtosecond laser and (b) picosecond
laser simulation. The peak laser intensity in pure $2\omega$ injecting scheme
is $I_{1}=0,I_{2}=65\%I_{0}$. The polarization direction of Laser 1 in
orthogonal and parallel polarization cases are $\hat{\bf e}_{1}=\hat{\bf
e}_{z}$ and $\hat{\bf e}_{1}=\hat{\bf e}_{y}$ separately, where the direction
of Laser 2 is kept at $\hat{\bf e}_{2}=\hat{\bf e}_{y}$.

According to Eq. (4), the strict harmonic relationship in dual-color effects
discussed in this article is not necessary, except that EIT mechanism requires
$\omega_{p}<\sqrt{\gamma}\omega_{2}<2\omega_{p},\omega_{2}-\omega_{1}<\omega_{p}/\sqrt{\gamma}$
Tie . We calculated the maximum of $\Delta a_{x}=a_{x,env}-a_{x,full}$ in Fig.
8(a), which is the difference between longitudinal fields excited by the full
field and laser envelope in Eq. (4). The frequency of Laser 2 varies from
$\omega_{2}=\omega_{1}/2$ to $\omega_{2}=\omega_{1}$. $(\Delta a_{x})_{max}$
is presented to evaluate the strength of longitudinal field modulation, and
the results show that the perturbation also takes place in anharmonic cases.
Figure 8(b) shows that a $10\%$ mixture of long-wavelength pulse in picosecond
injecting cases can lifting the high-energy part of spectrum significantly,
even though the frequency of Laser 2 is not an integral multiple of Laser 1.
The high tolerance of frequency relationship indicates that it could be
possible to obtain an extreme high-energy electron beam by injecting an
additional low-intensity but long-wavelength beam along with the main pulse.

Figure 8:  (a) The maximum of $\Delta a_{x}=a_{x,env}-a_{x,full}$ within
$k_{1}\xi\in[-50,0]$ calculated from Eq. (4) where $2a_{1}=a_{2}=5$ in Eq.
(5). (b) Energy spectrum of different injecting schemes, where
$I_{1}=I_{0},I_{2}=10\%I_{0}$ and $\lambda_{1}=1$ µm are fixed. The red and
green lines denote $\lambda_{2}=1.2$ µm ($\omega_{2}/\omega_{1}=5/6$) and
$\lambda_{2}=1.5$ µm ($\omega_{2}/\omega_{1}=2/3$) separately.

##  VI Conclusion

In summary, we have shown that the reservation of fundamental frequency
component can evidently lift the high-energy part of electron spectrum and
enhance the laser-plasma coupling process with different pulse duration (from
femtosecond to picosecond). The dual-color mixing effect are investigated and
attributed to the modulation of longitudinal field and the formation of
microchannel where different physical mechanisms such as relativistic EIT
would emerge with the second beam’s injection. The polarization direction and
frequency relationship required are also discussed. The results in this paper
can be applied to fast ignition schemes (for example, DCI) and experimental
designs of super-hot electron generation with different pulse durations.

###### Acknowledgements.

This work was supported by the Strategic Priority Research Program of Chinese
Academy of Sciences (Grant Nos. XDA25050300, XDA25010300, XDA25010100), the
National Key R&D Program of China (Grant No. 2018YFA0404801), and National
Natural Science Foundation of China (Grant No. 11827807 and 11775302), and the
Fundamental Research Funds for the Central Universities, the Research Funds of
Renmin University of China (20XNLG01). Computational resources have been
partially provided by the Physical Laboratory of High Performance Computing at
Renmin University of China.

## References

  * Zylstra  _et al._ (2022) A. B. Zylstra, O. A. Hurricane, D. A. Callahan, _et al._ , Nature 601, 542 (2022).
  * Abu-Shawareb  _et al._ (2022) H. Abu-Shawareb, R. Acree, P. Adams, _et al._ , Physical Review Letters 129, 75001 (2022).
  * Colaïtis  _et al._ (2022) A. Colaïtis, D. P. Turnbull, I. V. Igumenschev, D. Edgell, R. C. Shah, O. M. Mannion, C. Stoeckl, D. Jacob-Perkins, A. Shvydky, R. Janezic, A. Kalb, D. Cao, C. J. Forrest, J. Kwiatkowski, S. Regan, W. Theobald, V. N. Goncharov, and D. H. Froula, Physical Review Letters 129, 95001 (2022).
  * Schlossberg  _et al._ (2021) D. J. Schlossberg, G. P. Grim, D. T. Casey, A. S. Moore, R. Nora, B. Bachmann, L. R. Benedetti, R. M. Bionta, M. J. Eckart, J. E. Field, D. N. Fittinghoff, M. Gatu Johnson, V. Geppert-Kleinrath, E. P. Hartouni, R. Hatarik, W. W. Hsing, L. C. Jarrott, S. F. Khan, J. D. Kilkenny, O. L. Landen, B. J. Macgowan, A. J. Mackinnon, K. D. Meaney, D. H. Munro, S. R. Nagel, A. Pak, P. K. Patel, B. K. Spears, P. L. Volegov, and C. V. Young, Physical Review Letters 127, 125001 (2021).
  * Tabak  _et al._ (1994) M. Tabak, J. Hammer, M. E. Glinsky, W. L. Kruer, S. C. Wilks, J. Woodworth, E. M. Campbell, M. D. Perry, and R. J. Mason, Physics of Plasmas 1, 1626 (1994).
  * Kodama  _et al._ (2001) R. Kodama, P. A. Norreys, K. Mima, A. E. Dangor, R. G. Evans, H. Fujita, Y. Kitagawa, K. Krushelnick, T. Miyakoshi, N. Miyanaga, T. Norimatsu, S. J. Rose, T. Shozaki, K. Shigemori, A. Sunahara, M. Tampo, K. A. Tanaka, Y. Toyama, T. Yamanaka, and M. Zepf, Nature 412, 798 (2001).
  * Tabak  _et al._ (2006) M. Tabak, D. Hinkel, S. Atzeni, E. M. Campbell, and K. Tanaka, Fusion Science and Technology 49, 254 (2006).
  * Kitagawa  _et al._ (2022) Y. Kitagawa, Y. Mori, K. Ishii, R. Hanayama, S. Okihara, Y. Arikawa, Y. Abe, E. Miura, T. Ozaki, O. Komeda, H. Suto, Y. Umetani, A. Sunahra, T. Johzaki, H. Sakagami, A. Iwamoto, Y. Sentoku, N. Nakajima, S. Sakata, K. Matsuo, R. S. Mirfayzi, J. Kawanaka, S. Fujiokua, K. Tsubakimoto, K. Shigemori, K. Yamanoi, A. Yogo, A. Nakao, M. Asano, H. Shiraga, T. Motohiro, T. Hioki, and H. Azuma, Nuclear Fusion 62 (2022), 10.1088/1741-4326/ac7966.
  * Nagatomo  _et al._ (2019) H. Nagatomo, T. Johzaki, T. Asahina, M. Hata, Y. Sentoku, K. Mima, and H. Sakagami, Nuclear Fusion 59 (2019), 10.1088/1741-4326/ab3a5d.
  * Zhang  _et al._ (2020) J. Zhang, W. M. Wang, X. H. Yang, D. Wu, Y. Y. Ma, J. L. Jiao, Z. Zhang, F. Y. Wu, X. H. Yuan, Y. T. Li, and J. Q. Zhu, Philos. T. R. A 378, 20200015 (2020).
  * Rusby  _et al._ (2021) D. R. Rusby, P. M. King, A. Pak, N. Lemos, S. Kerr, G. Cochran, I. Pagano, A. Hannasch, H. Quevedo, M. Spinks, M. Donovan, A. Link, A. Kemp, S. C. Wilks, G. J. Williams, M. J. Manuel, Z. Gavin, A. Haid, F. Albert, M. Aufderheide, H. Chen, C. W. Siders, A. Macphee, and A. Mackinnon, Physical Review E 103, 1 (2021).
  * Robinson  _et al._ (2014) A. P. Robinson, D. J. Strozzi, J. R. Davies, L. Gremillet, J. J. Honrubia, T. Johzaki, R. J. Kingham, M. Sherlock, and A. A. Solodov, Nuclear Fusion 54 (2014), 10.1088/0029-5515/54/5/054003, arXiv:1304.1040 .
  * Nagatomo  _et al._ (2021) H. Nagatomo, T. Johzaki, M. Hata, Y. Sentoku, S. Fujioka, K. Mima, and H. Sakagami, Nuclear Fusion 61 (2021), 10.1088/1741-4326/ac2d5a.
  * Sawada  _et al._ (2020) H. Sawada, C. M. Salinas, F. N. Beg, H. Chen, A. J. Link, H. S. McLean, P. K. Patel, Y. Ping, and G. J. Williams, Plasma Physics and Controlled Fusion 62 (2020), 10.1088/1361-6587/ab8413.
  * Woldegeorgis  _et al._ (2019) A. Woldegeorgis, S. Herzer, M. Almassarani, S. Marathapalli, and A. Gopal, Physical Review E 100, 1 (2019).
  * Gibbon and Bell (1992) P. Gibbon and A. R. Bell, Physical Review Letters 68, 1535 (1992).
  * Forslund  _et al._ (1977) D. W. Forslund, J. M. Kindel, and K. Lee, Physical Review Letters 39, 284 (1977).
  * Wilks  _et al._ (1992) S. C. Wilks, W. L. Kruer, M. Tabak, and A. B. Langdon, Physical Review Letters 69, 1383 (1992).
  * Ivanov  _et al._ (2014) K. A. Ivanov, S. A. Shulyapov, P. A. Ksenofontov, I. N. Tsymbalov, R. V. Volkov, A. B. Savel’Ev, A. V. Brantov, V. Y. Bychenkov, A. A. Turinge, A. M. Lapik, A. V. Rusakov, R. M. Djilkibaev, and V. G. Nedorezov, Physics of Plasmas 21 (2014), 10.1063/1.4896348.
  * Yabuuchi  _et al._ (2010) T. Yabuuchi, B. S. Paradkar, M. S. Wei, J. A. King, F. N. Beg, R. B. Stephens, N. Nakanii, M. Hatakeyama, H. Habara, K. Mima, K. A. Tanaka, and J. T. Larsen, Physics of Plasmas 17 (2010), 10.1063/1.3447878.
  * Wang and Sheng (2008) W.-M. Wang and Z.-M. Sheng, Physics of Plasmas 15, 013101 (2008).
  * Faure  _et al._ (2006) J. Faure, C. Rechatin, A. Norlin, A. Lifschitz, Y. Glinec, and V. Malka, Nature 444, 737 (2006).
  * Dodd  _et al._ (2004) E. S. Dodd, J. K. Kim, and D. Umstadter, Physical Review E 70, 056410 (2004).
  * Chitgar  _et al._ (2020) Z. M. Chitgar, P. Gibbon, J. Böker, A. Lehrach, and M. Büscher, Physics of Plasmas 27 (2020), 10.1063/1.5117503, arXiv:1911.12762 .
  * Zeng  _et al._ (2015) M. Zeng, M. Chen, L. L. Yu, W. B. Mori, Z. M. Sheng, B. Hidding, D. A. Jaroszynski, and J. Zhang, Physical Review Letters 114, 1 (2015), arXiv:1406.7106 .
  * Yu  _et al._ (2014) L.-L. Yu, E. Esarey, C. B. Schroeder, J.-L. Vay, C. Benedetti, C. G. R. Geddes, M. Chen, and W. P. Leemans, Physical Review Letters 112, 125001 (2014).
  * Morace  _et al._ (2019) A. Morace, N. Iwata, Y. Sentoku, K. Mima, Y. Arikawa, A. Yogo, A. Andreev, S. Tosaki, X. Vaisseau, Y. Abe, S. Kojima, S. Sakata, M. Hata, S. Lee, K. Matsuo, N. Kamitsukasa, T. Norimatsu, J. Kawanaka, S. Tokita, N. Miyanaga, H. Shiraga, Y. Sakawa, M. Nakai, H. Nishimura, H. Azechi, S. Fujioka, and R. Kodama, Nature Communications 10, 1 (2019).
  * Zhou  _et al._ (2022) G. Zhou, W.-m. Wang, Y. Li, and J. Zhang, Physics of Plasmas 29, 052704 (2022).
  * Scott  _et al._ (2013) R. H. Scott, F. Pérez, M. J. Streeter, E. L. Clark, J. R. Davies, H. P. Schlenvoigt, J. J. Santos, S. Hulin, K. L. Lancaster, F. Dorchies, C. Fourment, B. Vauzour, A. A. Soloviev, S. D. Baton, S. J. Rose, and P. A. Norreys, New Journal of Physics 15 (2013), 10.1088/1367-2630/15/9/093021.
  * Baton  _et al._ (2008) S. D. Baton, M. Koenig, J. Fuchs, A. Benuzzi-Mounaix, P. Guillou, B. Loupias, T. Vinci, L. Gremillet, C. Rousseaux, M. Drouin, E. Lefebvre, F. Dorchies, C. Fourment, J. J. Santos, D. Batani, A. Morace, R. Redaelli, M. Nakatsutsumi, R. Kodama, A. Nishida, N. Ozaki, T. Norimatsu, Y. Aglitskiy, S. Atzeni, and A. Schiavi, Physics of Plasmas 15, 1 (2008).
  * Wang  _et al._ (2015) W. M. Wang, P. Gibbon, Z. M. Sheng, and Y. T. Li, Physical Review E - Statistical, Nonlinear, and Soft Matter Physics 91, 1 (2015), arXiv:1409.1808 .
  * Tajima and Dawson (1979) T. Tajima and J. M. Dawson, Physical Review Letters 43, 267 (1979).
  * Brantov  _et al._ (2008) A. V. Brantov, T. Z. Esirkepov, M. Kando, H. Kotaki, V. Y. Bychenkov, and S. V. Bulanov, Physics of Plasmas 15 (2008), 10.1063/1.2956989.
  * Gibbon (2005) P. Gibbon,  _Short Pulse Laser Interactions with Matter_, 1st ed. (Imperial College Press, London, 2005) p. 75.
  * Dalla and Lontano (1993) S. Dalla and M. Lontano, Physics Letters A 173, 456 (1993).
  * V.  _et al._ (1989) B. S. V., K. V. I., and S. A. S., JETP Letters 50, 176 (1989).
  * Sheng  _et al._ (2002) Z.-M. Sheng, K. Mima, Y. Sentoku, M. S. Jovanović, T. Taguchi, J. Zhang, and J. Meyer-ter Vehn, Physical Review Letters 88, 055004 (2002).
  * Meyer-ter Vehn and Sheng (1999) J. Meyer-ter Vehn and Z. M. Sheng, Physics of Plasmas 6, 641 (1999).
  * Huang  _et al._ (2013) T. W. Huang, C. T. Zhou, and X. T. He, Physical Review E 87, 053103 (2013).
  * Wan  _et al._ (2018) Y. Wan, C.-H. Pai, C. J. Zhang, F. Li, Y. P. Wu, J. F. Hua, W. Lu, C. Joshi, W. B. Mori, and V. Malka, Physical Review E 98, 013202 (2018).
  * Sprangle  _et al._ (1992) P. Sprangle, E. Esarey, J. Krall, and G. Joyce, Physical Review Letters 69, 2200 (1992).
  * Pukhov and Meyer-Ter-Vehn (1996) A. Pukhov and J. Meyer-Ter-Vehn, Physical Review Letters 76, 3975 (1996).
  * Pukhov  _et al._ (1999) A. Pukhov, Z.-M. Sheng, and J. Meyer-ter Vehn, Physics of Plasmas 6, 2847 (1999).
  * Arefiev  _et al._ (2012) A. V. Arefiev, B. N. Breizman, M. Schollmeier, and V. N. Khudik, Physical Review Letters 108, 1 (2012).
  * Kaw and Dawson (1970) P. Kaw and J. Dawson, Physics of Fluids 13, 472 (1970).
  * (46) Tie-Huai Zhang, Wei-Min Wang, Yu-Tong Li, and Jie Zhang, Electromagnetically Induced Transparency in Strongly Relativistic Regime, submitted .
  * Harris (1996) S. E. Harris, Physical Review Letters 77, 5357 (1996).
  * Gordon  _et al._ (2000a) D. F. Gordon, W. B. Mori, and C. Joshi, Physics of Plasmas 7, 3145 (2000a).
  * Gordon  _et al._ (2000b) D. F. Gordon, W. B. Mori, and C. Joshi, Physics of Plasmas 7, 3156 (2000b).
  * Pukhov and Meyer-Ter-Vehn (1997) A. Pukhov and J. Meyer-Ter-Vehn, Physical Review Letters 79, 2686 (1997).
  * Kemp  _et al._ (2008) A. J. Kemp, Y. Sentoku, and M. Tabak, Physical Review Letters 101, 8 (2008).
  * Weng  _et al._ (2012) S. M. Weng, M. Murakami, P. Mulser, and Z. M. Sheng, New Journal of Physics 14 (2012), 10.1088/1367-2630/14/6/063026.