Jerk and injection current implementations

doc/tutorial/tutorial.tex

@@ -2689,12 +2689,12 @@ \subsection*{Further ideas}
 
 \end{enumerate}
 
-\sectiontitle{25: Gallium Arsenide -- Nonlinear shift current}
+\sectiontitle{25: Gallium Arsenide -- Shift- and jerk-currents}
 
 \begin{itemize}
 
-\item Outline: \textit{Calculate the nonlinear shift current of inversion asymmetric fcc Gallium Arsenide. In preparation for this example it may be useful to read Ref.
-\cite{ibanez-azpiroz_ab_2018} }
+\item Outline: \textit{Calculate the nonlinear shift- and jerk-currents of inversion asymmetric fcc Gallium Arsenide. In preparation for this example it may be useful to read Ref.
+\cite{ibanez-azpiroz_ab_2018}, \cite{ruiz-ArXiv2023} }
 
 
 
@@ -2789,7 +2789,7 @@ \subsection*{Shift current $\sigma^{abc}$}
 
 Finally, \verb|sc_phase_conv| controls the phase convention used for the Bloch sums. 
 \verb|sc_phase_conv=1| uses the so-called tight-binding convention, whereby the Wannier centres are included
-into the phase, while   \verb|sc_phase_conv=2| leaves the Wannier centres out of the phase.
+into the phase, while  \verb|sc_phase_conv=2| leaves the Wannier centres out of the phase.
 These two possible conventions are explained in Ref. \cite{pythtb}. 
 Note that the overall shift-current spectrum does not depend on the chosen convention,
 but the individual terms that compose it do.
@@ -2807,6 +2807,42 @@ \subsection*{Shift current $\sigma^{abc}$}
 gnuplot> plot 'GaAs-sc_xyz.dat' u 1:2 w l
 \end{verbatim} 
 
+\subsection*{Jerk current $\iota^{abcd}$}
+The jerk-current tensor of GaAs has only three independent components that are finite, namely $\iota^{xxxx}, \iota^{xyyx}, \iota^{xxyy}$. For its computation, set
+\begin{verbatim}
+berry = true
+berry_task = jc
+\end{verbatim}
+Like the shift-current, the jerk-current is a frequency-dependent quantity. 
+The frequency window and step is controlled by \verb|kubo_freq_min|, \verb|kubo_freq_max| and 
+\verb|kubo_freq_step|, as explained in the users guide. 
+
+The jerk-current requires an integral over the Brillouin zone. The interpolated k-mesh is controlled by \verb|berry_kmesh|,
+which has been set to 
+\begin{verbatim}
+berry_kmesh = 100 100 100
+ \end{verbatim}
+We also need to input the value of the Fermi level in eV:
+\begin{verbatim}
+fermi_energy = [insert your value here]
+\end{verbatim}
+
+Due to possible degeneracies involved in the calculation of the second band derivatives required to compute the jerk-current tensor components, one needs to consider a small broadening parameter to avoid numerical problems (see parameter $\eta$ in Eq. (1) of Ref. \cite{ruiz-ArXiv2023} and related discussion).
+As in the case of the shift current, this parameter is controlled by \verb|sc_eta|. It is normally found that values between 0.01 eV and 0.1 eV 
+yield an stable spectrum. The default value is set to  $0.04$ eV.
+
+On output, the program generates a set of 36 files named \verb|SEED-jc_****.dat|, 
+which correspond to the different tensor components of the jerk-current
+(note that the 45 remaining components until totaling $3\times3\times3\times3=81$ 
+can be obtained from the 36 outputed by taking into account that $\iota^{abcd}$ is 
+symmetric under $a\leftrightarrow d$ and $b\leftrightarrow c$ index exchanges).
+For plotting the finite jerk-current components of GaAs $\iota^{xxxx}, \iota^{xyyx}, \iota^{xyxy}$ (units of Am/V$^{3}$s$^{2}$),
+\begin{verbatim}
+myshell> gnuplot
+gnuplot> plot 'GaAs-jc_xxxx.dat' u 1:2 w l
+gnuplot> replot 'GaAs-jc_xyyx.dat' u 1:2 w l
+gnuplot> replot 'GaAs-jc_xxyy.dat' u 1:2 w l
+\end{verbatim} 
 
 \sectiontitle{26: Gallium Arsenide -- Selective localization and constrained centres}
 
@@ -3795,6 +3831,100 @@ \subsection*{Expansion coefficients}
 
 \end{itemize}
 
+\sectiontitle{34: TaAs -- Injection-current tensor components}
+
+\begin{itemize}
+
+\item Outline: \textit{Calculate the nonlinear injection-current of TaAs. In preparation for this example it may be useful to read Ref. \cite{ruiz-ArXiv2023} }
+
+
+
+\item Directory: \verb|examples/example34/|
+
+\item Input files:
+
+\begin{itemize}
+
+\item[--] \verb|TaAs.scf| \textit{The {\tt PWSCF} input file for ground state calculation}
+\item[--] \verb|TaAs.nscf| \textit{The {\tt PWSCF} input file to obtain Bloch states on a uniform grid}
+\item[--] \verb|TaAs.pw2wan| \textit{The input file for} \verb|pw2wannier90|
+\item[--] \verb|TaAs.win| \textit{The} \verb|wannier90| \textit{and} \verb|postw90| \textit{input file}
+
+
+\end{itemize}
+
+
+\begin{enumerate}
+
+\item Run {\tt PWSCF} to obtain the ground state of Tantalum Arsenide
+
+\verb|pw.x < TaAs.scf > scf.out|
+
+
+\item Run {\tt PWSCF} to obtain the ground state of Tantalum Arsenide
+
+\verb|pw.x < TaAs.nscf > nscf.out|
+
+\item Run {\tt Wannier90} to generate a list of the required overlaps (written into the \verb|TaAs.nnkp| file)
+
+\verb|wannier90.x -pp TaAs|
+
+
+\item Run {\tt pw2wannier90} to compute:
+
+\begin{itemize}
+\item[--] The overlaps $\langle u_{n\bf{k}}|u_{n\bf{k+b}}\rangle$ between spinor 
+Bloch states (written in the \verb|TaAs.mmn| file) 
+\item[--] The projections for the starting guess (written in the \verb|TaAs.amn| file)
+
+\end{itemize}
+
+
+\verb|pw2wannier90.x < TaAs.pw2wan > pw2wan.out|
+
+\item Run {\tt wannier90} to  compute MLWFs
+
+\verb|wannier90.x TaAs|
+
+\item Run {\tt postw90} to compute expansion coefficients
+
+\verb|postw90.x TaAs| 
+
+
+\end{enumerate}
+
+\end{itemize}
+
+\subsection*{Injection current $\eta^{abc} = \eta^{abc}_S + \eta^{abc}_A$}
+The injection-current tensor of TaAs has only one independent component that is finite, namely $\eta^{xzx}_A$. For its computation, set
+\begin{verbatim}
+berry = true
+berry_task = ic
+\end{verbatim}
+Like the shift- or the jerk-current, the injection-current is a frequency-dependent quantity. 
+The frequency window and step is controlled by \verb|kubo_freq_min|, \verb|kubo_freq_max| and 
+\verb|kubo_freq_step|, as explained in the users guide. 
+
+The injection-current requires an integral over the Brillouin zone. The interpolated k-mesh is controlled by \verb|berry_kmesh|,
+which has been set to 
+\begin{verbatim}
+berry_kmesh = 100 100 100
+ \end{verbatim}
+We also need to input the value of the Fermi level in eV:
+\begin{verbatim}
+fermi_energy = [insert your value here]
+\end{verbatim}
+
+On output, the program generates a set of $18+9$ files named \verb|SEED-ic_S_***.dat| and \verb|SEED-ic_A_***.dat|, 
+which correspond to the different tensor components of the injection-current, which have been separated into symmetric ($S$) ($18$) and antisymmetric ($A$) ($9$)
+parts with respect to the $b\leftrightarrow c$ index exchange. The symmetric part is completely real, can only be nonzero in materials exhibiting time-reversal symmetry beaking 
+and only provides current in presence of linearly polarized light.
+The antisymmetric part is completely imaginary and only provides current in presence of circularly polarized light.
+For plotting imaginary part of the finite injection-current component of TaAs $\eta^{xzx}_A$ (units of A/V$^{2}$s),
+\begin{verbatim}
+myshell> gnuplot
+gnuplot> plot 'TaAs-ic_xzx.dat' u 1:2 w l
+\end{verbatim} 
 
 %\cleardoublepage
 

doc/user_guide/berry.tex

@@ -389,6 +389,85 @@ \section{{\tt berry\_task=sc}: shift current}
 
 Please cite Ref.~\cite{ibanez-azpiroz_ab_2018} when publishing shift-current results using this method.
 
+\section{{\tt berry\_task=ic}: injection current}
+
+The injection-current contribution to the second-order response
+is characterized by a frequency-dependent third-rank tensor \cite{ruiz-ArXiv2023}
+\begin{equation}\label{eq:injcurrent}
+\begin{split}
+\eta^{abc}(0;\omega,-\omega)=&-\frac{\pi e^3}{\hbar^2 \Omega_c N_k}
+\sum_{\bm{k}} \sum_{n,m}(f_{n\bm{k}}-f_{m\bm{k}})
+\times
+\Delta^a_{nm}(\bm{k})\left(r^b_{ nm}(\bm{k})r^c_{mn}(\bm{k})\right)\\
+&\times \delta(\omega_{mn\bm{k}}-\omega),
+\end{split}
+\end{equation}
+where $a,b,c$ are spatial indexes
+and $\omega_{mn\bm{k}}=(\epsilon_{n\bm{k}}-\epsilon_{m\bm{k}})/\hbar$.
+The expression in Eq.~\ref{eq:injcurrent} involves 
+the dipole matrix element, defined in Eq. \eqref{eq:r}
+and the difference between the band velocities
+\begin{equation}
+	\label{eq:delta}
+	\Delta^a_{nm} = \frac{\partial}{\partial k^a} \omega_{nm\bm{k}},
+\end{equation}
+which is calculated by using the scheme proposed in Ref. \cite{yates-prb07}.
+
+As in the case of the optical conductivity and the shift-current, a broadened delta function can be 
+applied in Eq.~\ref{eq:injcurrent} by means of the parameter
+$\eta$ (see Eq.~\ref{eq:lorentzian}) using the keyword 
+{\tt[kubo\_]smr\_fixed\_en\_width}, and adaptive smearing can 
+be employed using the keyword {\tt [kubo\_]adpt\_smr}. 
+
+The output files related to the injection-current take the form \verb|SEED-ic_S_***.dat| and \verb|SEED-ic_A_***.dat|. 
+The \verb|S| and \verb|A| labels denote the symmetric and antisymmetric parts of Eq. \eqref{eq:injcurrent} with respect to the $b\leftrightarrow c$ exchange, defined as
+\begin{subequations}
+	\label{eq:inj_sym_asym_components}
+	\begin{align}
+	\eta^{abc}_S(\omega) & = \frac{1}{2}\left[\eta^{abc}(\omega) + \eta^{acb}(\omega)\right],\label{eq:inj_sym}\\
+	\eta^{abc}_A(\omega) & = \frac{1}{2}\left[\eta^{abc}(\omega) - \eta^{acb}(\omega)\right].\label{eq:inj_asym}
+	\end{align}
+\end{subequations}
+The symmetric part is completely real and only nonvanishing in time-reversal lacking systems and the antisymmetric part is purely imaginary, as argued in \cite{ruiz-ArXiv2023}.
+
+Please cite Ref. \cite{ruiz-ArXiv2023} when publishing injection-current results using this method.
+
+\section{{\tt berry\_task=jc}: jerk current}
+
+The jerk-current contribution to the third-order response
+is characterized by a frequency-dependent fourth-rank tensor \cite{ruiz-ArXiv2023}
+\begin{equation}\label{eq:jerkcurrent}
+\begin{split}
+\iota^{abcd}(0;\omega,-\omega, 0)=&\frac{2\pi e^4}{\hbar^3 \Omega_c N_k}
+\sum_{\bm{k}} \sum_{n,m}(f_{n\bm{k}}-f_{m\bm{k}})
+\times
+\vartheta^{ad}_{nm}(\bm{k})\left(r^b_{ nm}(\bm{k})r^c_{mn}(\bm{k})\right)\\
+&\times \delta(\omega_{mn\bm{k}}-\omega),
+\end{split}
+\end{equation}
+where $a,b,c, d$ are spatial indexes
+and $\omega_{mn\bm{k}}=(\epsilon_{n\bm{k}}-\epsilon_{m\bm{k}})/\hbar$.
+The expression in Eq.~\ref{eq:injcurrent} involves 
+the dipole matrix element, defined in Eq. \eqref{eq:r}
+and the difference between the inverse effective masses
+\begin{equation}
+	\label{eq:eff_mass}
+	\vartheta^{ad}_{nm} = \frac{\partial^2}{\partial k^a \partial k^d} \omega_{nm\bm{k}},
+\end{equation}
+which is calculated by using the scheme proposed in Ref. \cite{yates-prb07}.
+The jerk-current photoconductivity tensor $\iota^{abcd}$ in \eqref{eq:jerkcurrent} is symmetric under the $b\leftrightarrow c$ and $a\leftrightarrow d$ exchanges.
+
+As in the case of the optical conductivity, the shift-current and injection-current, a broadened delta function can be 
+applied in Eq.~\ref{eq:injcurrent} by means of the parameter
+$\eta$ (see Eq.~\ref{eq:lorentzian}) using the keyword 
+{\tt[kubo\_]smr\_fixed\_en\_width}, and adaptive smearing can 
+be employed using the keyword {\tt [kubo\_]adpt\_smr}. 
+
+The output files related to the jerk-current take the form \verb|SEED-jc_****.dat|. 
+
+Please cite Ref. \cite{ruiz-ArXiv2023} when publishing jerk-current results using this method.
+
+
 \section{{\tt berry\_task=kdotp}: $k\cdot p$ coefficients}
 \label{sec:kdotp}
 

doc/wannier90.bib

@@ -743,3 +743,12 @@ @article{Lihm_shift_eta_2021
   url = {https://link.aps.org/doi/10.1103/PhysRevB.103.247101}
 }
 
+@article{ruiz-ArXiv2023,
+author = {R. Puente{-}Uriona, {\'A}lvaro and S. Tsirkin, Stepan and Souza, Ivo and Iba\~nez{-}Azpiroz, Julen},
+title = "{\textit{Ab initio} study of the nonlinear optical properties and d.c. photocurrent of the Weyl semimetal TaIrTe$_4$}",
+journal = {ArXiv e-prints},
+eprint = {2302.03090},
+year = {2023},
+url = {https://doi.org/10.48550/arXiv.2302.03090},
+}
+