Critical Laser Intensity of Phase-Matched High-Order Harmonic Generation in Noble Gases

Abstract

The efficient generation of high-order harmonic radiation has been a challenging task since the early days of strong-field physics. An essential requirement to achieve efficient high-order harmonic generation inside a gas medium is the phase matching of the high-order harmonic radiation and the incident laser pulse. The dominant contribution to the wave–vector mismatch $\Delta k$ is associated with the ionization probability of the medium. In this work, we derive two analytical formulas to calculate the critical intensity of a general linearly polarized laser pulse that obey the phase-matching condition $\Delta k=0$. The analytic formulas are valid in the tunneling regime (ADK model) and the regime of the tunnel and multi-photon ionization (PPT model), respectively. We compare our results to numerical computations and discuss the scaling of the critical intensity depending on the pulse duration and the wavelength of a realistic incident laser pulse. The analytical approach demonstrated in this work is highly accurate and can compete with the existing numerical computational methods by an error of less than 1% and a decrease in the computation time of approximately 4 to 6 orders of magnitude. This enables complex theoretical studies of the efficiency scaling in HHG or to consider the effects of ground state depletion efficiently.

1. Introduction

High-order harmonic generation (HHG) is an essential technique to produce coherent ultra-short light pulses from tabletop-sized devices. The radiation from these devices can be used in industry and in fundamental scientific research, which highlights them as a promising source for a broad range of applications.

The macroscopic conversion efficiency of the high-order harmonic (HH) process is on the order of ${10}^{-4}$ at $30\mathrm{eV}$ [1] down to ${10}^{-9}$ at $300\mathrm{eV}$ [2]. Specifically, the conversion efficiency of HHG from near- to mid-infrared laser pulses decreases exponentially with increasing wavelengths, which has been demonstrated in experimental and theoretical work [3,4,5]. This low efficiency induced the interest of the strong-field community to find specific experimental conditions that optimize the HH yield [6,7,8,9,10]. The highest experimental conversion efficiencies can be achieved by simultaneously optimizing macroscopic effects such as absorption [11], laser focusing [12,13], and phase matching [14,15,16].

Phase matching (PM) is one of the essential requirements to efficiently produce HH radiation. Therefore, the phase velocity of the incident laser pulse and the generated HH radiation need to match while propagating through the gas medium. The PM of complex pulse structures, e.g., twisted light [16,17], is often relatively complicated; thus, linearly polarized laser pulses are advantageous to simplify the PM. Commonly used media for the generation of HH radiation are isotropic noble gases, which are also the focus of this work. The conventional approach of second or third-order harmonic generation of using birefringence to phase-match the fundamental and harmonic beam cannot be used. However, it can be shown that the contribution of plasma (free electrons + ions) and neutral gas dispersion, the geometrical phase, and the intrinsic phase can be balanced to achieve a transient phase-matching window [14].

In this work, we analytically investigate the PM conditions for a partially ionized noble gas in the free-focusing regime that can be generalized to other geometries (e.g., tight focusing or waveguide). We explicitly show that the critical peak intensity of a realistic laser pulse depends inverse logarithmically on the pulse duration and, further, has a nontrivial wavelength dependency associated with the refraction index. In the discussed experimental setup, other PM contributions are nullified by each other. That enables us to discuss an experimental setup under realistic conditions and demonstrate a newly developed calculation technique. However, this calculation technique is not limited to a specific experimental setup and may be adapted, e.g., in complex analytical or numerical studies, where a numerically calculated ionization rate is not suited.

This paper is structured as follows: In Section 2, we derive the analytical formulas of the critical field intensities for which PM is achieved at the peak of the incident laser pulse. In more detail, in Section 2.1, we revisit the standard model of PM for HHG and specify the PM condition of this work. In Section 2.2, we derive and discuss the formula of the critical intensity in the tunneling regime. Section 2.3 generalizes the derived formula to the intermediate regime where the tunnel and multi-photon ionization needs to be considered, which incorporates many realistic scenarios. In Section 3, we compare our results with numerical calculations and discuss the parameter range for which the developed method is valid. Finally, we conclude our findings in Section 4. In the following, we use atomic units ($\hslash =e=m_e=4\pi {ϵ}_0=1$) unless stated otherwise.

2. Model and Method

In Figure 1, we can see a standard HHG setup with the incident laser pulse on the left (orange-red) and the emitted HH radiation on the right (blue). This work focuses on the linear propagation effects of the HH radiation and the incident laser pulse inside the gas medium. In particular, we focus on the phase-matching of these light fields, which is highlighted in Figure 1a by the dashed circle.

The medium consists of an arbitrary noble gas characterized by its ionization potential $I_p$ and orbital angular momentum quantum number ℓ of the respective atomic state. We assume that the incident laser pulse is a paraxial linearly polarized plane wave propagating in the z-direction (optical axis) that irradiates the gas target (free-focusing regime). The optical axis is perpendicular to the target surface in the x-y plane while the target is centered one Rayleigh range behind the focus of the incident laser pulse $z=z_R$ (not explicitly shown in Figure 1a). Experimentally, this is done to favor the phase-matched generation of the short trajectories [18]. In the free-focusing regime, the variation of the electric field strength on the optical axis can be neglected for gas targets that are much thinner than the Rayleigh range such that the atoms in the interaction region are affected by a constant field strength. In the following Section 2.1, we use the SI unit system.

2.1. Phase Matching of High-Order Harmonic Radiation

For a coherent build-up of the HH radiation, the wave–vector mismatch has to be minimized, resulting in efficient HHG. The wave–vector mismatch is defined as

$$\begin{array}{c}\hfill \Delta k=k\left(q\omega \right)-qk\left(\omega \right),\end{array}$$

(1)

where the first term is associated with the HH radiation and the second with the incident laser pulse. The wave–vector mismatch consists of four major contributions in the free-focusing regime, as discussed in [19]

$$\begin{array}{c}\hfill \Delta k=\Delta k_{atom}+\Delta k_{plasma}+\Delta k_{Gouy}+\Delta k_{intrinsic}\equiv 0,\end{array}$$

(2)

where $\Delta k_{atom}$ and $\Delta k_{plasma}$ consider the mismatch associated with the neutral and ionized atoms, respectively. The term $\Delta k_{Gouy}$, on the other hand, denotes a geometric contribution by the Gouy phase, and $\Delta k_{intrinsic}$ describes the wave–vector mismatch related to the intrinsic phase of the emitted photons. The Gouy wave–vector mismatch has a positive contribution, whereas the intrinsic phase results in a negative contribution to the total wave–vector mismatch. These two contributions nullify; eventually $\Delta k_{Gouy}+\Delta k_{intrinsic}\approx 0$, as the target is centered one Rayleigh length behind the focus of the laser pulse [18,20,21]. This leaves us with the atomic and plasma contributions that can be combined to an dispersion mismatch as

$$\begin{array}{cc}\hfill & \Delta k_{disp}=\Delta k_{atom}+\Delta k_{plasma},\hfill \\ \hfill \Rightarrow & -q\frac{\omega }{c_0}\frac{\rho }{N_{atm}}\Delta \delta \left(1-\frac{\eta }{{\eta }_c}\right)=0,\hfill \end{array}$$

(3)

where $\omega$ is the frequency of the laser pulse, $c_0$ is the speed of light in the vacuum, $N_{atm}$ is the particle number at standard conditions, $\rho$ is the particle density, and $\Delta \delta =n_0\left(\omega \right)-n_0\left(q\omega \right)$ is the difference between the refractive indices of the laser pulse and the qth harmonic under standard conditions. The respective refraction indices are calculated as in [22,23,24]. The ionization probability $\eta$ denotes macroscopically the relative amount of atoms in an ionized state compared to the total number of atoms. Furthermore, ${\eta }_c$ is the so-called critical ionization probability that is defined as

$$\begin{array}{c}\hfill {\eta }_c={\left(1+\frac{2\pi r_e{N}_{atm}{c}_0^2}{{\omega }^2\Delta \delta }\right)}^{-1},\end{array}$$

(4)

where $r_e=e^2/\left(4\pi {ϵ}_0{m}_e{c}_0^2\right)$ is the classical electron radius. To fulfill (3), the term inside the parentheses should vanish, which is fulfilled if the ionization probability matches the critical ionization probability $\eta ={\eta }_c$ in Equation (3). The critical ionization probability (4) of commonly used noble gases is shown in Figure 1b for a fundamental laser wavelength of $\lambda =1000\mathrm{nm}$.

2.2. ADK Model: Critical Field Intensity

In this work, we investigate strong-field ionization in the context of phase-matched high-order harmonic generation. Apart from solving the time-dependent Schrödinger equation numerically, the YI model [25] does describe the strong-field ionization sufficiently well. As HHG is known to operate in the regime of tunnel and multi-photon ionization, we restrict our investigations to approximations of the YI model that incorporate the respective ionization processes (ADK and PPT).

This section is separated into two parts. In the first part, we sketch the idea of a solution that obeys Equation (3). The second part incorporates a detailed derivation as well as a discussion of the applied approximations.

2.2.1. Critical Field Intensity as Phase Matching Condition

For a pulsed laser, it is impossible to fulfill Equation (3) in general because the ionization probability increases monotonically over time. Therefore, the PM condition should be achieved at the peak of the laser pulse, resulting in the highest number of ionized electrons. The ADK ionization rate is defined as [20,26]

$$\begin{array}{cc}\hfill w^{ADK}\left(t\right)=& \frac{1}{(2\ell +1)}\sum _{m=-\ell }^{\ell }{w}_m^{ADK}\left(t\right),\hfill \\ \hfill =& \sum _{m=-\ell }^{\ell }\frac{|C_{n^{*}{l}^{*}}{|}^2{G}_{lm}{I}_p}{(2\ell +1)}{\left(\frac{2{\left(2I_p\right)}^{3/2}}{F\left(t\right)}\right)}^{2n^{*}-\left|m\right|-1}{e}^{-\frac{2{\left(2I_p\right)}^{3/2}}{3F\left(t\right)}},\hfill \\ \hfill \equiv & \sum _{m=-\ell }^{\ell }{\kappa }_m{\left(\frac{F_0}{F\left(t\right)}\right)}^{g_m+1}{e}^{-\frac{F_0}{F\left(t\right)}},\hfill \end{array}$$

(5)

with the atomic species dependent factor ${\kappa }_m=3^{g_m+1}{\left|C_{n^{*}{l}^{*}}\right|}^2{G}_{lm}{I}_p{(2\ell +1)}^{-1}$, which defines the time-independent amplitude of the ionization rate concerning the initial state of the ionized electron. The parameters $|C_{n^{*}{l}^{*}}{|}^2$ and $G_{lm}$ are defined in Equations (4.65) and (4.66) in [20], respectively. Moreover, we define the terms

$$\begin{array}{c}\hfill F_0=\frac{2}{3}{\left(2I_p\right)}^{\frac{3}{2}},\mathrm{and}{g}_m=2n^{*}-\left|m\right|-2,\end{array}$$

(6)

as atom-specific constants, and the electric field of the laser pulse reads

$$\begin{array}{c}\hfill F\left(t\right)=E_{0}f\left(t\right)cos\left(\omega t\right),\end{array}$$

(7)

$$\begin{array}{c}\hfill f\left(t\right)=e^{-4ln\left(2\right){\left(\frac{t}{\tau }\right)}^2}.\end{array}$$

(8)

Here, $f\left(t\right)$ denotes a Gaussian envelope with a full-width half maximum pulse duration of $\tau$. The ADK ionization probability in the peak of the laser pulse then reads [20]

$$\begin{array}{c}\hfill \eta ={\eta }^{ADK}=1-e^{-{\int }_{-\infty }^0{w}^{ADK}\left(t\right)dt}.\end{array}$$

(9)

One can see the dependency of the ionization probability on the laser frequency $\omega$, the atomic ionization potential $I_p$, the pulse duration $\tau$, and the peak intensity of the laser pulse $I=\sqrt{E_0}$. In the experimental setup, it is often not possible to adapt the atomic species or the wavelength of the laser pulse to optimize these parameters with regard to the PM condition (3). On the other hand, the intensity and the pulse duration can be controlled relatively easily compared to other parameters. Following this idea, we can define the critical intensity $I_c$ of the laser pulse as a function of the wavelength, atomic species, and pulse duration that fulfills (3) with

$$\begin{array}{c}\hfill \eta (I,\tau ,I_p,\lambda )={\eta }_c\Rightarrow \eta \left(I_c(\tau ,I_p,\lambda )\right)={\eta }_c.\end{array}$$

(10)

Numerical algorithms are often used to find the critical intensity $I_c$ from Equation (9) that fulfills the phase matching condition in Equation (3). These algorithms integrate the ADK ionization rate in Equation (5) (or PPT rate in (25)) for a given intensity and iteratively optimize them until the intensity converges to the critical intensity $I=I_c$. In other words, these algorithms have the temporal integration of Equation (9) inside a root-finding algorithm. This algorithm can be computationally heavy because the numerical integration needs to be executed for each iteration of the root-finding algorithm.

Within our analytical approach, which is in detail derived in Section 2.2.2, the critical intensity of the ADK model $I_c\to I_c^{\left(0\right)}$ is approximated by

$$\begin{array}{cc}\hfill I_c^{\left(0\right)}(\tau ,I_p,\lambda )& ={\left[\frac{g_0}{F_0}{W}_0\left(-\frac{1}{g_0}{\left(\frac{D_0}{\tau }\right)}^{1/g_0}\right)\right]}^{-2},\mathrm{for}{g}_0<0,\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill I_c^{\left(0\right)}(\tau ,I_p,\lambda )& ={\left[\frac{g_0}{F_0}{W}_{-1}\left(-\frac{1}{g_0}{\left(\frac{D_0}{\tau }\right)}^{1/g_0}\right)\right]}^{-2},\mathrm{for}{g}_0>0,\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill D_0& =2\sqrt{2ln\left(2\right)}\frac{|ln(1-{\eta }_c\left)\right|}{{\kappa }_0}.\hfill \end{array}$$

(13)

where ${\kappa }_m$ is defined below Equation (5) and $g_0$ in Equation (6) with $m\to 0$. The Lambert W function is denoted by $W_i$ for its respective real branch with $i\in \{0,-1\}$. For all neutral noble gases, except Xenon, the parameter $g_0$ is smaller than zero such that Equation (11) can be used for helium, neon, argon, and krypton, whereas Equation (12) can be used for xenon. Note that the variables $F_0$ and $g_0$ depend only on the atomic species, whereas $D_0\propto ln(1-{\eta }_c\left(\lambda \right))$ depends on the wavelength of the incident laser pulse. The dependency of the critical intensity on the laser pulse is therefore imprinted in the critical ionization probability ${\eta }_c\to {\eta }_c\left(\lambda \right)$ and the pulse duration $\tau$. The atomic-species dependent parameters for the noble gases can be found in

Table 1. List of atom-dependent parameters to calculate the critical intensity in Equations (11) and (12).

	He	Ne	Ar	Kr	Xe
$g_0$	−0.5122	−0.4114	−0.1417	−0.0283	0.1182
$F_0$	1.6196	1.3303	0.8311	0.6958	0.5612
${\kappa }_0$	6.5714	6.4204	6.1200	6.0226	5.9116

.

$$\begin{array}{cc}\hfill \frac{|ln(1-{\eta }_c\left({\lambda }_0\right)\left)\right|}{{\tau }_0}=\frac{|ln(1-{\eta }_c^0\left)\right|}{{\tau }_0}& \equiv \frac{|ln(1-{\eta }_c^1\left)\right|}{{\tau }_1}\hfill \\ \hfill \Rightarrow {\tau }_1={\tau }_0\frac{|ln(1-{\eta }_c^1\left)\right|}{|ln(1-{\eta }_c^0\left)\right|}& \approx {\tau }_0\frac{{\eta }_c^1}{{\eta }_c^0}\mathrm{for}{\eta }_c^0,{\eta }_c^1<<1\hfill \end{array}$$

(14)

Equation (14) denotes a condition for the pulse duration ${\tau }_1$ in which the critical intensity of the incident laser pulse remains invariant while changing the parameters from ${\lambda }_0,{\tau }_0\to {\lambda }_1,{\tau }_1$. The factor $|ln(1-{\eta }_c\left)\right|$ is physically equivalent to the time-integrated ionization rate at the peak of the laser pulse

$$\begin{array}{c}\hfill |ln(1-{\eta }_c)|={\int }_{\infty }^{0}dt{w}^{ADK}\left(t\right).\end{array}$$

(15)

If we reformulate Equation (14), we can find the following invariance relation:

The critical field intensity is invariant under a parameter change from $({\tau }_0,{\lambda }_0)$ to $({\tau }_1,{\lambda }_1)$, if the relative time integrated ionization rate equals the relative pulse duration.

$$\begin{array}{c}\hfill \frac{{\tau }_1}{{\tau }_0}=\frac{n_1{\lambda }_1}{n_0{\lambda }_0}=\frac{ln(1-{\eta }_c\left({\lambda }_1\right))}{ln(1-{\eta }_c\left({\lambda }_0\right))}=\frac{ln(1-{\eta }_c^1)}{ln(1-{\eta }_c^0)}\mathrm{and}\mathrm{if}{\eta }_c^1,{\eta }_c^0<<1\Rightarrow \frac{{\tau }_1}{{\tau }_0}\approx \frac{{\eta }_c^1}{{\eta }_c^0}\end{array}$$

(16)

Simplifying (11) and (12) even further reveals a simple scaling law of the critical intensity for both positive and negative $g_0$ [27].

$$\begin{array}{cc}\hfill I_c^{\left(0\right)}(\tau ,I_p,\lambda )& \approx {\left[\frac{g_0}{F_0}ln\left(\frac{1}{|g_0|}{\left(\frac{D_0}{\tau }\right)}^{1/|g_0|}\right)\right]}^{-2}\propto F_0^2{ln}^{-2}\left(\frac{D_0}{\tau }\right)\hfill \\ \hfill \Rightarrow I_c^{\left(0\right)}\left(\tau \right)& \propto {ln}^{-2}\left(\frac{|ln(1-{\eta }_c\left)\right|}{\tau }\right)\hfill \end{array}$$

(17)

With this fairly simple scaling law, we are able to describe the general proportionality of the critical intensity for all noble gases with respect to the pulse duration and the critical ionization probability of a monochromatic Gaussian laser pulse.

To summarize, the critical intensity calculated in (11) depends on the incident laser pulse via the pulse duration $\tau$ and, indirectly, on the wavelength of the laser pulse ${\eta }_c\left(\lambda \right)$, where it obeys the proportionality $I_c^{\left(0\right)}\propto {ln}^{-2}\left(\frac{|ln(1-{\eta }_c\left)\right|}{\tau }\right)$. The wavelength dependency is directly related to the critical ionization probability in which the refractive index is sensitive to the wavelength of the incident laser pulse. The critical intensity does not significantly depend on electrons ionized from states with the magnetic quantum number $m=\pm 1$. At the same time, the dependency on the ionization potential $I_p$ cannot be derived so easily.

2.2.2. Derivation of the Critical Intensity

In the following derivation, we make use of a number of approximations. All of these approximations can be reduced to one of the following major assumptions.

• $F_0/E_0>>1$: The ionization potential $F_0\propto I_p^{3/2}$ is much larger than the peak amplitude of the incident laser pulse $E_0$.
• $n⪆10$: The number of optical cycles of the laser pulse is sufficiently large.
• $m\ne \pm 1$: Contributions to the ionization rate from atomic states with magnetic quantum number $m=\pm 1$ can be neglected.

To derive the formula of $I_c\equiv I_c(\tau ,I_p,\lambda )$ we focus on the temporal integration in (5) and pull the summation over the magnetic quantum number m out of the integral.

$$\begin{array}{c}\hfill \frac{1}{2\ell +1}{\int }_{-\infty }^{0}dt{w}_m^{ADK}\left(t\right)=\frac{1}{2(2\ell +1)}{\int }_{-\infty }^{\infty }{w}_m^{ADK}=\frac{{\kappa }_m}{2}{\int }_{-\infty }^{\infty }dt{\left(\frac{F_0}{F\left(t\right)}\right)}^{g_m+1}{e}^{-\frac{F_0}{F\left(t\right)}}\end{array}$$

(18)

In the next step, we approximate the Gaussian envelope by assuming a parametric dependency on the time such that the amplitude of the envelope is constant within each half cycle $t\in \left[-\frac{T_0}{4},\frac{T_0}{4}\right]$. This approximation is valid for a slowly varying envelope, which implies a sufficiently large number of optical cycles.

$$\begin{array}{cc}\hfill \frac{1}{2\ell +1}{\int }_{-\infty }^{0}dt{w}_m^{ADK}\left(t\right)\approx & \frac{{\kappa }_m}{2}\sum _{k=-\infty }^{\infty }{\int }_{-T_0/4}^{T_0/4}dt{\left(\frac{F_0}{E_0{f}_{k}cos\left(\omega t\right)}\right)}^{g_m+1}{e}^{-\frac{F_0}{E_0{f}_{k}cos\left(\omega t\right)}}\hfill \\ \hfill \approx & {\kappa }_m{\left(\frac{F_0}{E_0}\right)}^{g_m+1}\sum _{k=-\infty }^{\infty }{\int }_0^{T_0/4}dt{e}^{-\frac{F_0}{E_0{f}_{k}cos\left(\omega t\right)}}\hfill \end{array}$$

(19)

In the last line, we neglected the factor $f_{k}cos\left(\omega t\right)$ in the denominator, as the dominant contribution is associated with the argument of the exponential function. To integrate the remaining expression, we further approximate the integral as

$$\begin{array}{cc}\hfill {\int }_0^{T_0/4}dt{e}^{-\frac{F_0}{E_0{f}_{k}cos\left(\omega t\right)}}=& \frac{T_0}{2\pi }{\int }_0^{\pi /2}dt{e}^{-\frac{F_0}{E_0{f}_{k}cos\left(x\right)}}\approx \frac{T_0}{2\pi }{\int }_0^{\infty }dt{e}^{-\frac{F_0}{E_0{f}_k}cosh\left(x\right)},\hfill \\ \hfill \Rightarrow {\int }_0^{T_0/4}dt{e}^{-\frac{F_0}{E_0{f}_{k}cos\left(\omega t\right)}}\approx & \frac{T_0}{2\pi }{K}_0\left(\frac{F_0}{E_0{f}_k}\right),\hfill \end{array}$$

(20)

with the modified Bessel function of order zero $K_0\left(\frac{F_0}{E_0{f}_k}\right)$. Using the series expansion of the modified Bessel function for large $x=\frac{F_0}{E_0{f}_k}$ ($K_0\left(x\right)\approx \sqrt{\frac{\pi }{2x}}{e}^{-x}$), we obtain

$$\begin{array}{cc}\hfill \frac{1}{2\ell +1}{\int }_{-\infty }^{0}dt{w}_m^{ADK}\left(t\right)\approx & \frac{T_0{\kappa }_m}{2\pi }{\left(\frac{F_0}{E_0}\right)}^{g_m+1}\sqrt{\frac{\pi }{2}}{\left(\frac{F_0}{E_0}\right)}^{-1/2}\sum _{k=-\infty }^{\infty }{e}^{-\frac{F_0}{E_0{f}_k}}.\hfill \end{array}$$

(21)

Note that the contributions of the envelope in the denominator are neglected. The remaining summation can be executed as

$$\begin{array}{cc}\hfill \sum _{k=-\infty }^{\infty }{e}^{-\frac{F_0}{E_0{f}_k}}=& \sum _{k=-\infty }^{\infty }{e}^{-\frac{F_0}{E_0}{e}^{(ln\left(2\right){\left(k/n\right)}^2)}}\approx \sum _{k=-\infty }^{\infty }{e}^{-\frac{F_0}{E_0}(1+ln\left(2\right){\left(k/n\right)}^2)},\hfill \\ \hfill \approx & e^{-\frac{F_0}{E_0}}{\int }_{-\infty }^{\infty }dk{e}^{-\frac{F_0}{E_0}ln\left(2\right){\left(k/n\right)}^2}=n\sqrt{\frac{\pi }{ln\left(2\right)}}{\left(\frac{F_0}{E_0}\right)}^{-1/2}{e}^{-\frac{F_0}{E_0}},\hfill \\ \hfill \Rightarrow \sum _{k=-\infty }^{\infty }{e}^{-\frac{F_0}{E_0{f}_k}}\approx & n\sqrt{\frac{\pi }{ln\left(2\right)}}{\left(\frac{F_0}{E_0}\right)}^{-1/2}{e}^{-\frac{F_0}{E_0}}.\hfill \end{array}$$

(22)

Here, we used the definition of the pulse duration as $\tau =n{T}_0$ for an integer number n of optical cycles. Inserting (22) into (21) yields us an analytic formula of the ionization rates with

$$\begin{array}{c}\hfill \frac{1}{2\ell +1}{\int }_{-\infty }^{0}dt{w}_m^{ADK}\left(t\right)\approx \frac{\tau {\kappa }_m}{2\sqrt{2ln\left(2\right)}}{\left(\frac{F_0}{E_0}\right)}^{g_m}{e}^{-\frac{F_0}{E_0}}.\end{array}$$

(23)

By inserting specific parameters, one can see that the contribution of the terms $w_{\pm 1}^{ADK}$ are low and can be neglected such that we can replace $\eta \to {\eta }_c$ and write (9) as

$$\begin{array}{cc}\hfill -ln(1-{\eta }_c)\approx & \frac{1}{(2\ell +1)}{\int }_{-\infty }^0{w}_0^{ADK}\left(t\right)dt,\hfill \\ \hfill {\left(\frac{F_0}{\sqrt{I_c}}\right)}^{g_0}{e}^{-\frac{F_0}{\sqrt{I_c}}}=& 2\sqrt{2ln\left(2\right)}\frac{|ln(1-{\eta }_c\left)\right|}{\tau {\kappa }_0}\equiv \frac{D_0}{\tau }.\hfill \end{array}$$

(24)

Solving this equation for the critical intensity $I_c$ yields us the solution of (3) with the ansatz (10) by

$$\begin{array}{cc}\hfill I_c^{\left(0\right)}(\tau ,I_p,\lambda )& ={\left[\frac{g_0}{F_0}{W}_0\left(-\frac{1}{g_0}{\left(\frac{D_0}{\tau }\right)}^{1/g_0}\right)\right]}^{-2},\mathrm{for}{g}_0<0,\hfill \\ \hfill I_c^{\left(0\right)}(\tau ,I_p,\lambda )& ={\left[\frac{g_0}{F_0}{W}_{-1}\left(-\frac{1}{g_0}{\left(\frac{D_0}{\tau }\right)}^{1/g_0}\right)\right]}^{-2},\mathrm{for}{g}_0>0.\hfill \\ \hfill D_0& =2\sqrt{2ln\left(2\right)}\frac{|ln(1-{\eta }_c\left)\right|}{{\kappa }_0}.\hfill \end{array}$$

2.3. PPT: Critical Field Intensity

The analytical method to calculate the critical intensity discussed in the last section can be extended to incorporate the mixed regime in which tunnel and multi-photon ionization need to be considered simultaneously. This regime is usually well described by the PPT model [20,28]. Unfortunately, the PPT ionization rates have a more complex dependency on the laser intensity than the ADK rates. To proceed further, we apply a perturbative approach. The applied technique is similar to the one demonstrated in Section 2.2 such that we refer to the specific equations if necessary.

The PPT ionization rate is defined as

$$\begin{array}{cc}\hfill w^{PPT}\left(t\right)=& \frac{1}{(2\ell +1)}\sum _{m=-\ell }^{\ell }{w}_m^{PPT}\left(t\right),\hfill \\ \hfill =& \sum _{m=-\ell }^{\ell }{\kappa }_m{\left(\frac{F_0}{F\left(t\right)}\right)}^{g_m+1}{e}^{-\frac{F_0}{F\left(t\right)}g\left(\gamma \right)}\hfill \\ & \hfill \times \frac{4}{\sqrt{3\pi }}\frac{1}{\left|m\right|!}\frac{{\gamma }^2}{1+{\gamma }^2}{\left(1+{\gamma }^2\right)}^{\left(\right|m|+1)/2}\sum _{q\ge q_{thr}}^{\infty }{A}_q(\omega ,\gamma ),\end{array}$$

(25)

with the unknown parameters defined in [20] next to Equation (4.64). To find an approximation of the critical intensity $I_c\to I_c^{\left(1\right)}$ concerning the PPT ionization rates, we assume that the Keldysh parameter can be approximated (if $\gamma \approx 1$) by the critical intensity of the ADK rates such that

$$\begin{array}{c}\hfill \gamma =\omega \sqrt{\frac{2I_p}{I_c}}\approx \gamma \left(I_c^{\left(0\right)}\right)\equiv {\gamma }^{\left(0\right)}.\end{array}$$

(26)

Inserting ${\gamma }^{\left(0\right)}$ into (25) can then be written as

$$\begin{array}{c}\hfill w^{PPT}\left(t\right)\approx \sum _{m=-\ell }^{\ell }{\kappa }_m^{\left(0\right)}{\left(\frac{F_0^{\left(0\right)}}{F\left(t\right)}\right)}^{g_m+1}{e}^{-\frac{F_0^{\left(0\right)}}{F\left(t\right)}},\end{array}$$

(27)

where we define

$$\begin{array}{cc}\hfill {\kappa }_m^{\left(0\right)}=& \frac{{\delta }_{l,\left|m\right|}+{\delta }_{\ell -1,\left|m\right|}}{\left|m\right|!}{\left(g\left({\gamma }^{\left(0\right)}\right)\right)}^{-(g_m+1)}\hfill \\ & \times {\kappa }_m\frac{4}{\sqrt{3\pi }}\frac{{\left({\gamma }^{\left(0\right)}\right)}^2}{1+{\left({\gamma }^{\left(0\right)}\right)}^2}{\left(1+{\left({\gamma }^{\left(0\right)}\right)}^2\right)}^{\left(\right|m|+1)/2}\sum _{q\ge q_{thr}}^{\infty }{A}_q(\omega ,{\gamma }^{\left(0\right)}),\hfill \end{array}$$

(28)

$$\begin{array}{cc}\hfill F_0^{\left(0\right)}=& F_{0}g\left({\gamma }^{\left(0\right)}\right).\hfill \end{array}$$

(29)

The Kronecker deltas ensure a correct summation over the magnetic quantum number m for $\ell \le 1$. The total ionization rate (25) then results in

$$\begin{array}{cc}\hfill w^{PPT}\left(t\right)=& \left[{\kappa }_0^{\left(0\right)}{\left(\frac{F_0^{\left(0\right)}}{F\left(t\right)}\right)}^{g_0+1}+2{\kappa }_1^{\left(0\right)}{\left(\frac{F_0^{\left(0\right)}}{F\left(t\right)}\right)}^{g_1+1}\right]e^{-\frac{F_0^{\left(0\right)}}{F\left(t\right)}},\hfill \\ \hfill =& \left[{\kappa }_0^{\left(0\right)}+2{\kappa }_1^{\left(0\right)}{\left(\frac{F_0^{\left(0\right)}}{F\left(t\right)}\right)}^{-1}\right]{\left(\frac{F_0^{\left(0\right)}}{F\left(t\right)}\right)}^{g_0+1}{e}^{-\frac{F_0^{\left(0\right)}}{F\left(t\right)}},\hfill \\ \hfill \approx & \left[{\kappa }_0^{\left(0\right)}+2{\kappa }_1^{\left(0\right)}{\left(\frac{F_0^{\left(0\right)}}{E_0^{\left(0\right)}}\right)}^{-1}\right]{\left(\frac{F_0^{\left(0\right)}}{F\left(t\right)}\right)}^{g_0+1}{e}^{-\frac{F_0^{\left(0\right)}}{F\left(t\right)}},\hfill \end{array}$$

(30)

$$\begin{array}{cc}\hfill {\kappa }^{\left(0\right)}=& \left[{\kappa }_0^{\left(0\right)}+2{\kappa }_1^{\left(0\right)}{\left(\frac{F_0^{\left(0\right)}}{E_0^{\left(0\right)}}\right)}^{-1}\right]\hfill \end{array}$$

(31)

$$\begin{array}{cc}\hfill \Rightarrow w^{PPT}\left(t\right)\approx & {\kappa }^{\left(0\right)}{\left(\frac{F_0^{\left(0\right)}}{F\left(t\right)}\right)}^{g_0+1}{e}^{-\frac{F_0^{\left(0\right)}}{F\left(t\right)}}.\hfill \end{array}$$

(32)

In the first line of the equation, we made use of the inversion symmetry of the magnetic quantum number in the PPT model such that $w_1^{PPT}\left(t\right)=w_{-1}^{PPT}\left(t\right)$. From the second to the third line, we assumed that the temporal contribution of the electric field $F\left(t\right)$ in the rectangular brackets is small. Furthermore, we assumed that the electric field strength can be approximated sufficiently well by the critical electric field strength $E_0^{\left(0\right)}=\sqrt{I_c^{\left(0\right)}}$ calculated from the ADK ionization rates. One can see that the mathematical structure of (32) is similar to the integrand in (18).

Using the calculation technique discussed in Section 2.2.2 yields the critical intensity calculated by the PPT ionization rates as

$$\begin{array}{cc}\hfill I_c^{\left(1\right)}(\tau ,I_p,\lambda )& ={\left[\frac{g_0}{F_0^{\left(0\right)}}{W}_0\left(-\frac{1}{g_0}{\left(\frac{D^{\left(0\right)}}{\tau }\right)}^{1/g_0}\right)\right]}^{-2},\mathrm{for}{g}_0<0,\hfill \end{array}$$

(33)

$$\begin{array}{cc}\hfill I_c^{\left(1\right)}(\tau ,I_p,\lambda )& ={\left[\frac{g_0}{F_0^{\left(0\right)}}{W}_{-1}\left(-\frac{1}{g_0}{\left(\frac{D^{\left(0\right)}}{\tau }\right)}^{1/g_0}\right)\right]}^{-2},\mathrm{for}{g}_0>0.\hfill \end{array}$$

(34)

$$\begin{array}{cc}\hfill D^{\left(0\right)}=& 2\sqrt{2ln\left(2\right)}\frac{|ln(1-{\eta }_c\left)\right|}{{\kappa }^{\left(0\right)}}.\hfill \end{array}$$

(35)

Note that we included the contributions for the magnetic quantum numbers $m=\pm 1$ such that we only used Assumptions 1 and 2 expressed at the beginning of Section 2.2.2 to analytically calculate the critical intensity within the PPT model. The scaling behavior of $I_c^{\left(1\right)}$ is more complex than the one of the critical intensity $I_c^{\left(0\right)}$ (11), as $F_0^{\left(0\right)}$ and ${\kappa }^{\left(0\right)}$ both depend on the critical intensity $I_c^{\left(0\right)}$ itself. Therefore, one may hardly see some general scaling behavior. Nevertheless, because we applied a perturbative approach, the general behavior should not differ strongly from the critical intensity $I_c^{\left(0\right)}$. This means that the dependency on the intensity in ${\kappa }^{\left(0\right)}$ and $F_0^{\left(0\right)}$ should be sufficiently small if the Keldysh parameter obeys $\gamma ⪅1$.

If the Keldysh parameter is expected to fulfill $\gamma ⪆1$, the method demonstrated in this section can be applied again to increase the accuracy.

$$\begin{array}{cc}\hfill \gamma \left(I_c^{\left(0\right)}\right)\equiv {\gamma }^{\left(0\right)}& \to \gamma \left(I_c^{\left(1\right)}\right)\equiv {\gamma }^{\left(1\right)}\hfill \\ \hfill \Rightarrow F_0^{\left(0\right)}\to F_0^{\left(1\right)};D^{\left(0\right)}& \to D^{\left(1\right)};I_c^{\left(1\right)}\to I_c^{\left(2\right)}\hfill \end{array}$$

(36)

In fact, this procedure can be applied iteratively such that the ith iteration delivers the $(i+1)$th critical intensity as

$$\begin{array}{cc}\hfill I_c^{(i+1)}(\tau ,I_p,\lambda )& ={\left[\frac{g_0}{F_0^{\left(i\right)}}{W}_0\left(-\frac{1}{g_0}{\left(\frac{D^{\left(i\right)}}{\tau }\right)}^{1/g_0}\right)\right]}^{-2},\mathrm{for}{g}_0<0,\hfill \end{array}$$

(37)

$$\begin{array}{cc}\hfill I_c^{(i+1)}(\tau ,I_p,\lambda )& ={\left[\frac{g_0}{F_0^{\left(i\right)}}{W}_{-1}\left(-\frac{1}{g_0}{\left(\frac{D^{\left(i\right)}}{\tau }\right)}^{1/g_0}\right)\right]}^{-2},\mathrm{for}{g}_0>0.\hfill \end{array}$$

(38)

$$\begin{array}{cc}\hfill D^{\left(i\right)}=& 2\sqrt{2ln\left(2\right)}\frac{|ln(1-{\eta }_c\left)\right|}{{\kappa }^{\left(i\right)}}.\hfill \end{array}$$

(39)

This procedure can be applied until the critical intensity has converged. Typically, a convergence is achieved for $i\in \left[2,8\right]$, which in detail depends on the system itself $(\tau ,I_p,\lambda )$.

3. Results and Discussion

In the following section, we discuss the accuracy of the derived analytical formulas concerning the critical intensity and the closely related Keldysh parameter. Moreover, we define the parameter space for which the formulas of the critical intensity are valid. To do so, we define three methods to measure the error of our formulas

$$\begin{array}{cc}\hfill \Delta I_c^{\left(0\right)}& =\frac{|I_c^{\left(0\right)}-I_c^{ADK}|}{I_c^{ADK}},\hfill \end{array}$$

(40)

$$\begin{array}{cc}\hfill \Delta I_c^{\left(i\right)}& =\frac{|I_c^{\left(i\right)}-I_c^{PPT}|}{I_c^{PPT}},\mathrm{for}i>0,\hfill \end{array}$$

(41)

$$\begin{array}{cc}\hfill \Delta {\gamma }^{\left(i\right)}& =\frac{|{\gamma }^{\left(i\right)}-{\gamma }^{PPT}|}{{\gamma }^{PPT}},\mathrm{for}i>0,\hfill \end{array}$$

(42)

where the first two errors indicate the modulus of the difference between the exact numerical solution of the critical intensity ($I_c^{ADK}$ and $I_c^{PPT}$) and their respective analytic approximations $I_c^{\left(0\right)}$ and $I_c^{\left(i\right)}$. The error of the Keldysh parameter is defined to demonstrate the convergence of Equation (37). All calculations are performed for an HH photon energy of $25\mathrm{eV}$ unless stated otherwise.

In the following discussion, we refer to the critical intensities of the ADK and PPT models by (11) and (37), respectively. This restriction denotes a shorthand notation and does not exclude the case for $g_0>0$. Furthermore, different atomic species can be visually distinguished by different line colors and different incident wavelengths by a respective marker or line type.

3.1. Accuracy Critical Intensity: Tunnel Ionization (ADK)

We begin with a visual comparison of the critical intensity calculated analytically and numerically. The critical intensity in Figure 2a is shown about the number of optical periods n, which is proportional to the pulse duration $n=\tau /T_0$. Moreover, the wavelength of the incident laser pulse is fixed as $\lambda =1000\mathrm{nm}$ and the atomic species is denoted by the line color. The analytical results of Equation (11) are denoted by squares and the numerical results are as a solid line.

Because the difference between the numerical and analytical results is hardly visible, we display the error (40) of the critical intensity $I_c^{\left(0\right)}$ for a set of wavelengths $\lambda =1000,2000,3000\mathrm{nm}$ (line type) and noble gases (line color) in Figure 2b. It is quite remarkable that the error is below $1\%$ for all parameter configurations, especially if we consider the number of approximations and assumptions we made in our derivation.

The error is not significantly sensitive to the harmonic energy; thus, this error estimate is valid for HH photon energies in the range from $25\mathrm{eV}$ to $200\mathrm{eV}$. The upper limit of our calculations is arbitrarily set to $200\mathrm{eV}$, which can be seen in Figure 1b. Nevertheless, the critical ionization probability remains fairly constant for the laser pulse parameters discussed in this paper such that the upper limit can be easily extended to photon energies of $2\mathrm{keV}$.

One can see that the error of the critical intensity decreases with an increasing wavelength of the incident laser pulse as well as exponentially decreases with an increasing pulse duration. In addition, note that the error decreases with increasing ionization potential for all noble gases, except for helium. Because helium in the ground state does not have electrons in the p-orbital, the contributions of ionized electrons from bound states with $m=\pm 1$ do not contribute to the numerical calculations compared to the other noble gases.

The accuracy of Equation (11) compared to a numerical computation can be seen for various configurations in

Table 2. Accuracy of the critical intensity in the ADK model (40). The table is showing the lower limit of the parameters ${\lambda }_{\min },n_{\min },{\tau }_{\min }\propto {\lambda }_{\min }{n}_{\min }$ (wavelength, number of optical cycles, pulse duration) for which the error of the critical intensity is smaller than $\Delta I_c^{\left(0\right)}$. The error decreases for any configuration of higher parameters. We fixed the HH photon energy to $50\mathrm{eV}$ such that the error of the critical intensity remains roughly constant for higher photon energies and changes reasonably slowly for lower photon energies, such as the critical ionization probability in Figure 1b.

Element	${\mathit{\lambda }}_{\mathbf{min}}$ in nm	${\mathit{n}}_{\mathbf{min}}$	${\mathit{\tau }}_{\mathbf{min}}$ in fs	$\mathit{\Delta }{\mathit{I}}_{\mathit{c}}^{\left(0\right)}$
He	800	7	18.7	$\le 1.0\%$
	1600	10	53.4	$\le 0.7\%$
	3200	5	53.4	$\le 0.6\%$
Ne	800	3	8.0	$<0.01\%$
	1600	3	16.0	$<0.01\%$
	3200	3	32.0	$<0.01\%$
Ar	800	3	8.0	$<0.5\%$
	1600	3	16.0	$<0.4\%$
	3200	3	32.0	$<0.3\%$
Kr	800	4	10.7	$<0.7\%$
	1600	4	21.4	$<0.5\%$
	3200	4	42.7	$<0.4\%$
Xe	800	5	13.3	$\le 1.0\%$
	1600	5	26.7	$\le 0.7\%$
	3200	5	53.4	$\le 0.5\%$

. Here, we show the shortest wavelength ${\lambda }_{\min }$ and the smallest number of optical periods $n_{\min }$ for which the error is smaller than $\Delta I_c^{\left(0\right)}$. This holds for $n>n_{\min }$ and $\lambda >{\lambda }_{\min }$, as the error decreases for increasing n and $\lambda$. As a brief conclusion, the error of Equation (11) obeys

$$\begin{array}{c}\hfill \tau \ge 20\mathrm{fs},\lambda \ge 800\mathrm{nm},n\ge 8\Rightarrow \Delta I_c^{\left(0\right)}\le 1.0\%,\end{array}$$

(43)

and holds true for all noble gases. To compare the computation times, we fix the parameters $\lambda =1000\mathrm{nm},I_p=15.76\mathrm{eV}$ (Ar: $\ell =1$), and ${\eta }_c=0.02$ for a high-harmonic photon energy of $E_q>100\mathrm{eV}$. We compare the critical intensity for pulse duration of $n=20\tau /T_0$ and $n=100\tau /T_0$. The numerical computation is executed in the Julia programming language [29] and uses the QuadGK package [30] for numerical integration and the NLsolve package [31] for the root-finding algorithm. The resulting computation times are denoted as $t_A\left(n\right)$ and $t_N\left(n\right)$ for the analytical and numerical computation, respectively. The average computation times for $n=20\tau /T_0$ are $t_A\left(20\right)=736\mathrm{ns}$ compared to $t_N\left(20\right)=830·{10}^5\mathrm{ns}$; thus, the analytical method is roughly ${10}^5$ times faster. For the pulse duration $n=100\tau /T_0$, the analytical method is roughly ${10}^6$ times faster with $t_A\left(100\right)=725\mathrm{ns}$ compared to $t_N\left(100\right)=447·{10}^6\mathrm{ns}$. Here, one can see the potential relevance for more advanced investigations that need to calculate the critical intensity many times. Overall, the ratio for the computation times of the critical intensities, calculated from the respective PPT model, is ${10}^4$ to ${10}^5$. This ratio is approximately one order of magnitude lower than the ratio for the critical intensity in the ADK model, which may be explained by the iterative application of the critical intensity Equation (37).

For wavelengths shorter than $800\mathrm{nm}$, the ADK model is no longer valid. We proceed further with the critical intensity calculated within the PPT model.

3.2. Accuracy Critical Intensity: Tunnel and Multi-Photon Ionization (PPT)

Going on to shorter wavelengths of the laser pulse leads us to the regime in which tunnel and multi-photon ionization need to be considered simultaneously. In this regime, the Keldysh parameter is on the order of one such that we need to calculate the critical intensities with the PPT ionization model. We use Equation (37), where i increases with increasing Keldysh parameter until the Keldysh parameter converges.

In Figure 3a, one can see the error of the critical intensity (37) for $i=2$, similar to Figure 2b. In contrast to Figure 2b, the error of heavy noble gases such as krypton and xenon deviates from the exponential decrease with increasing pulse duration. In addition, the error gets fairly large for long pulse durations. This unusual scaling of the error in krypton and xenon is related to the number of iteration steps i in (37). Because the critical Keldysh parameter ${\gamma }_c^{\left(i\right)}$ is larger for small ionization potentials, as in the case of krypton and xenon, the number of iterations needs to be larger as well. Figure 3b shows the convergence of the Keldysh parameter with iteration step i on the x-axis and the error of the Keldysh parameter on the y-axis. Here, we fixed the pulse duration to $\tau =200T_0$. The iteration step $i=0$ calculates the Keldysh parameter by using the critical intensity of the ADK model. The Keldysh parameter converges for increasing iteration steps i and wavelengths of $\lambda =500,750,1000\mathrm{nm}$ (marker type) fairly quickly to an error of less than $1\%$.

Similar to the critical intensity in the ADK model, the error is not significantly sensitive to the energy of the HH photon, hence we can neglect it. The accuracy of the critical intensity in the PPT model is shown for various configurations in

Table 3. Accuracy of the critical intensity in the PPT model (41). The table is showing the lower limit of the parameters with the respective wavelength ${\lambda }_{\min },$ number of optical cycles $n_{\min }$, and the resulting pulse duration ${\tau }_{\min }\propto {\lambda }_{\min }{n}_{\min }$ for which the error of the critical intensity is smaller than $\Delta I_c^{\left(0\right)}$. The iteration parameter i is not fixed and is chosen such that the critical intensity is converged. We fixed the HH photon energy to $50\mathrm{eV}$ such that the critical intensity remains roughly constant for higher photon energies and, thus, the change in the error for lower photon energies is reasonably small.

Element	${\mathit{\lambda }}_{\mathbf{min}}$ in nm	${\mathit{n}}_{\mathbf{min}}$	$\mathit{\tau }$ in fs	$\mathit{\Delta }{\mathit{I}}_{\mathit{c}}^{\left(\mathit{i}\right)}$
He	250	5	4.2	$\le 2.5\%$
	515	5	8.6	$<1.5\%$
	800	10	26.7	$\le 1.0\%$
Ne	250	9	7.5	$\le 2.5\%$
	515	5	8.6	$\le 1.5\%$
	800	8	21.4	$\le 1.0\%$
Ar	250	10	8.3	$\le 2.5\%$
	515	3	5.2	$\le 1.5\%$
	800	4	10.7	$\le 1.0\%$
Kr	250	5	4.2	$<2.5\%$
	515	3	5.2	$<1.5\%$
	800	3	8.0	$<1.0\%$
Xe	250	3	2.5	$<2.5\%$
	515	3	5.2	$<1.0\%$
	800	3	8.0	$<1.0\%$

. Overall the accuracy of (37) obeys

$$\begin{array}{cc}\hfill \tau \ge 8.3\mathrm{fs},\lambda \ge 250\mathrm{nm},n\ge 10& \Rightarrow \Delta I_c^{\left(i\right)}\le 2.5\%,\hfill \end{array}$$

(44)

$$\begin{array}{cc}\hfill \tau \ge 8.6\mathrm{fs},\lambda \ge 515\mathrm{nm},n\ge 5& \Rightarrow \Delta I_c^{\left(i\right)}\le 1.5\%,\hfill \end{array}$$

(45)

$$\begin{array}{cc}\hfill \tau \ge 26.7\mathrm{fs},\lambda \ge 800\mathrm{nm},n\ge 10& \Rightarrow \Delta I_c^{\left(i\right)}\le 1.0\%,\hfill \end{array}$$

(46)

for all noble gases. Here, the computation time of the numerical integration is also 3 to 4 orders of magnitude higher than for our analytical approach (37). Note, the iteration step parameter i is chosen such that the critical intensity is converged and is therefore not fixed.

3.3. Critical Intensity: Comparison ADK—PPT

As we have already discussed the accuracy of the critical intensities (11) and (37), in this section, we focus on the interpretation of their respective results. In particular, we focus on the general scaling of the critical intensity concerning the pulse duration and the incident wavelength. Figure 4 shows two density plots of the critical intensity $I_c^{\left(i\right)}(\lambda ,\tau )$ (color scale) in which the variation of the wavelength is located on the horizontal axis and the variation of the pulse duration on the vertical axis. Figure 4a,b show the critical intensity of the ADK and the PPT model in argon, respectively. Specifically, the critical intensity within the ADK model increases constantly with shrinking wavelength and pulse duration. This is highlighted by the contour lines (solid) that denote parameter compositions ($\lambda ,\tau$) with constant critical intensity. These lines represent the solutions of Equation (16) inserted into $I_c^{\left(i\right)}(\lambda ,\tau )$. The pulse duration scaling matches well with our analytic approximation in Equation (17), whereas the wavelength scaling is related to the critical ionization probability that does not show a closed form of its scaling. As noted before, Figure 4b shows the critical intensity in the PPT model. First, let us focus on the dashed inset on the right side. Here, the scaling of the critical intensity matches well with the ADK model such that the respective structure of the contour lines agrees. However, outside of the inset, the scaling changes drastically. The deviation is remarkable, as the global scaling changes rapidly. The change in the global scaling is visualized by the purple stars to the left of the dashed inset. The position of these stars can be interpreted as follows.

If we select an arbitrary pair of parameters (${\lambda }_0,{\tau }_0$) inside the dashed inset, then we can constantly decrease the wavelength while keeping the pulse duration constant, which leads to a monotonic increase in the critical intensity. At a specific wavelength, the critical intensity will not increase anymore but remains constant and will decrease afterward. The purple star denotes the wavelength ${\lambda }_0$ at which the intensity remains constant such that it describes the transition between two separated regimes. In the first regime, the critical intensity increases for decreasing wavelength; in the second, the critical intensity decreases with decreasing wavelength. Even though the number of purple stars that are shown in Figure 4b is finite, in general, they represent an infinite collection of points.

From a more pragmatic point of view, the purple stars denote the highest critical intensity that can be applied to generate phase-matched high-order harmonics, for a given pulse duration $\tau$.

As noted above, the data shown in Figure 4b was computed by the analytical formula of the critical intensity in the PPT model Equation (37). The computation time for this $1000\times 1000$ data grid (wavelength and the pulse duration) is, on average, $97\mathrm{s}$. A numerical algorithm as defined below Equation (10) would need the same computation, approximately $97\mathrm{s}·{10}^4=270\mathrm{h}$, or 11 days.

4. Conclusions

In this work, we demonstrated an approach to calculate the critical intensity of phase-matched high-order harmonic generation analytically. The derived formulas are valid in the tunnel ionization regime (ADK-model) and get extended to the intermediate regime of tunnel and multi-photon ionization (PPT model). The analytic formulas for the critical intensity in both ionization regimes are highly accurate and can replace numerical computations while benefiting from a decrease in the computation time of approximately 4 to 6 orders of magnitude. We analyzed the accuracy of the developed approximations for the commonly used parameter space of all noble gases and listed the errors for specific sub-spaces. In addition, we discussed the analytic scaling properties of the critical intensity with regard to the pulse duration and wavelength in the tunnel-ionization regime. Finally, we highlighted the differences between the critical intensities in the respective ionization regimes and displayed the suppression of this intensity in the intermediate ionization regime.

Overall, we demonstrated that the developed formulas of the critical intensity are highly accurate and show a significantly decreased computation time compared to a numerical solution.