A New Way to Predict the Efficiency of Optical Limiters without Providing an Experiment: Processing of TDDFT Calculations for the Case of Pure Two-Photon Absorption.

Laser-power-limiting devices play a predominant role in photonics because of their potential for protecting human eyes and optical devices that are sensitive to intense laser beams. This paper describes a new methodology for predicting the efficiency of optical limiting based on electric-field-induced changes in absorption spectra calculated with the TDDFT quantum-chemical method. Analytical equations are derived to evaluate the optical thresholds and speed of switching on, the dynamic range, and the degree of nonlinear attenuation of the radiation fluxes for the case of two-photon absorption. Thus, the researcher does not need to conduct costly experiments to evaluate the suitability of nonlinear absorbers for the creation of optical limiters. The possibility of developing a forecasting model is demonstrated by an example of a series of stable slipped-cofacial phthalocyanine J-type dimers, which were synthesized and investigated previously.

Introduction

The optical limiting effect is due to the nonlinear attenuation of the power, energy, or radiation flux transmitted by an optical system to a certain fixed value regardless of the magnitude of the input signal.[1,2] At the same time, the absorbing material of the optical limiter must have high transmittance at low input radiation powers.[3,4] This result is commonly used in the development of devices for protecting sensitive optical sensors from damage by laser radiation and controlling or stabilizing laser radiation in various optical logical systems.[5−7] Optical limiting is caused by nonlinear effects in the interaction of intense light with an absorbing medium, such as nonlinear absorption, refraction, scattering, and optically induced phase transitions.[8,9] Among these mechanisms, nonlinear absorption is considered when organic materials with an extended π-system (for example, phthalocyanines and related dyes) are considered since it directly affects the transmission, causing it to decrease at high levels of input radiation intensity.[10] This is because the laser radiation leads to saturation of excited states, and the dye molecule exhibits the so-called reverse saturable absorption (RSA) effect.[11,12] For phthalocyanines, this effect is observed in the range of 400–600 nm, in which most known lasers operate. The best optical limiter should have a low power-on threshold and attenuate radiation to a fixed level.

In the experiment, the optical limiting effect is observed as a decrease in transmission on the growth of the radiation power, when the sample is placed at the focus of a collecting lens.[13,14] At the same time, it was shown that the optical effect observed in laser experiments can be reproduced using ab initio calculations.[15] For this purpose, TDDFT calculations are performed by perturbing the ground state with electric fields polarized in all Cartesian directions.[16] Thus, the electrical component of laser radiation is associated with a series of electric fields in simulation. The calculated absorption spectra of isolated organic dyes under different electric fields demonstrate significant variation in the spectral region of transparency. The simulated optical response as the “transmittance vs. field strength” relationship shows a tendency toward optical limiting similar to what is usually observed in practice.[15−17] To bring theory and practice closer together, primarily, one should compare the oscillation strength of the TDDFT transitions and the absorption coefficient of the dyes in solution or a matrix. Then, one can operate with such macroscopic parameters as linear transmission and optical layer thickness. Finally, to calculate the parameters responsible for the efficiency of optical limiting, it is necessary to derive a series of analytical mathematical expressions. In the future, this will make it possible to utilize chemical structures to predict the efficiency of laser radiation limiting without conducting costly experiments.

This article is devoted to solving the above-mentioned theoretical problem for the case of pure two-photon absorption (2PA), which we previously implemented to calculate the nonlinear absorption coefficient and optical limiting response for the series of stable slipped-cofacial phthalocyanine J-type dimers (Figure ).[18]

Figure 1

Common chemical structure of slipped-cofacial phthalocyanine J-dimers.

Methods

Quantum-Chemical Calculations

Initially, structures of the dyes 1a–f,[16,18] in which the tert-butyl substituents were replaced with hydrogen atoms, were fully optimized using density functional theory (DFT) implemented in the quantum-chemical program PRIRODA[19] without solvent effects and symmetry constraints. The gradient-corrected exchange-correlation Perdew, Burke, and Ernzerhof (PBE) functional,[20] as well as the cc-pVDZ basis set[21] were utilized for this purpose. The efficient resolution of identity and parallel implementation of evaluating both Coulomb and exchange-correlation integrals with optimized fitting Gaussian basis sets in the PRIRODA code permits the performance of calculations of the molecular systems with a large number of basis functions.[19,21] A large integration grid (which comprises about 8 × 105 points over calculated molecules) with a 5 × 10–8 accuracy parameter of the adaptively generated grid was used. This parameter is responsible for the precision of the exchange-correlation energy per atom. The 10–6 threshold on the orbital gradients at the energy calculations tag and 10–5 threshold on the molecular gradient at the geometry optimization procedure were employed. The valence shells were described by basis sets with the following contraction schemes: {6s2p}/[2s1p] on H; {10s7p3d}/[3s2p1d] on C, N, and O; {14s11p6d}/[5s4p2d] on Mg; and {19s15p11d5f}/[6s5p3d1f] on Co, Ni, Cu, and Zn atoms, respectively. The vibrational analysis confirmed whether the optimized geometries correspond to local minima without imaginary frequencies.

Further, the level of theory for the optimized dimers 1a–f was raised by including Grimme’s dispersion correction. Pure DFT GGA functional B97-D and 6-31 + G(d) basis sets were utilized to perform calculations on the GAMESS-US Software[22] to re-optimize the target structures. Functional B97-D belongs to one of the most accurate general-purpose GGAs. The performance for noncovalently bound systems including many pure van der Waals complexes is exceptionally good, reaching on the average CCSD(T) accuracy.[23] In turn, the 6-31 + G(d) basis set contains diffuse and polarization functions that are required to determine the convergence of the finite field property.[24] The polarizable continuum model (PCM) was used to control the solvent effect computations (solvent – THF). The optimized structures and general computed data are presented in Section .

To study the influence of the electric fields on the electronic transitions, the linear-response time-dependent density functional theory (LR-TDDFT) implemented in GAMESS-US Software was utilized. Thus, each molecule is placed in a constant in time and uniform in space external electric field, which is included in the Schrödinger equation for the ground state,[25] followed by the calculation of the energies of excited states and oscillator strengths. Here, we assume that the action of electric fields on molecules is short-lived in time and does not cause geometric changes. The described methodology subsequently allowed us to obtain strict correlations of theoretical and experimental data regarding the parameters of optical limiting.

The resulted TDDFT spectra for each absolute value of field were obtained by the simple accumulation of all transitions in an array, followed by executing a fit with the Lorentz function L(x) according to the formulain which w is the line half-width, and x is the position of the electronic transitions. The variable line half-width (in the range of 0.05–0.3 eV), stretch coefficient (1.15–1.2), and a displacement factor (+ca. 25 nm) were used for fitting the calculated spectra to the experimental ones[18] by applying the Spectroscopy module for EasyQuanto.[26] This is necessary to most accurately predict the shape of the spectra under the action of electric fields of different intensities on the dye molecules. All quantum chemical calculations were performed on an Intel/Linux cluster (Joint Supercomputer Centre of the Russian Academy of Sciences – http://jscc.ru). Visualization of the optimized structures (Figure ) was performed with the Mercury program obtained from the Cambridge Crystallographic Data Centre (http://ccdc.cam.ac.uk).

Figure 4

DFT (B97-D/6-31 + G(d))-optimized structures of phthalocyanine J-type dimers. tert-Butyl substituents were replaced with hydrogen atoms to reduce calculation time.

Analytical Calculations

Conversion of Electric Field to Power

The absolute value of the field strength in the quantum chemical programs is given in a.u. (1 a.u. = 5.14 × 1011 V m–1). The series of fields F was associated with the electric component of model laser radiation and converted to power units P (Wcm–2) according to the following formulawhere c is the speed of light in a vacuum, n is the refractive index of the medium (for THF, n = 1.408), ε0 is the vacuum permittivity (ε0 = 8.85 × 10–12 F m–1), F is the field intensity (in V m–1), dA is the average diameter of the aggregate absorbing the light, and dL is the laser spot diameter that is focused on the sample. Thus, we imagine an aggregation model when a homogeneous electrical induction encompasses a sphere in which the absorber molecules are enclosed and under the assumption that the aggregate exhibits the properties of a single isolated molecule. Further, we will assume that the energy of the model radiation is concentrated in a focused laser beam that passes through the aggregate, modifying its spectral characteristics. This hypothesis has allowed correlating the simulation results with the experiment.

Processing TDDFT Spectra

The intensities of the electronic transitions in the framework of TDDFT are expressed as oscillator forces, fosc, theoretically related to the probability of these transitions occurring. To compare the fosc values with intensities of the experimental bands, we introduced a macroscopic coefficient DF, the density factor of the medium, according to the following equationin which T0 is the linear transmission, α is the linear absorption coefficient (cm–1),

and d is the layer thickness (cm).

The density factor is a dimensionless quantity. For each model structure, DF was computed by an iterative procedure. First, a sequence of DF values from 1 to 0 was formed with a step of 0.02. Then, the absorption curves (eq ) were calculated to evaluate α in the framework of the pure 2PA absorption model (see Section ). In this way, a correlation between DF and T0 was established. The density factor was assumed to be constant for all sets of electric fields in the series.

For the numerical calculations, analytical wavelength and T0 were set to 532 nm and 0.7, respectively, as in the experiment.[18] In this way, we can establish the relationship between the calculated TDDFT spectra and the molar concentration of the sample.

Theoretical Background

Since the relationship between absorption and the external fields can be nonlinear, we used the following nonlinear functionin which a, b, and c are the fitting parameters.

This is a modified Cauchy’s distribution[27] represented as the smooth, continuous, and infinitely differentiable function for all x > 0. This function allows deriving all the required analytical expressions in a simple and understandable form. However, for the 2PA absorption model, a linear relationship between absorption and radiation power is assumed.[18] Therefore, we expand k(x) into the Taylor series, truncating high-order terms:where β is the effective nonlinear two-photon absorption coefficient.[28,29] The proposed approach allows one to derive equations for three-photon (3PA), mixed (2 + 3) PA, and other absorption models of any complexity. The optical limiting properties of dyes 1a–f were previously investigated in the framework of the 2PA absorption model at the 532 nm wavelength by applying nanosecond laser radiation (16 ns pulses).[18] In this way, we can establish the correlation between simulation and experiment (Section .).

The optical limiting curves (output signal) relate the output radiation flux with the input radiation, according to the equationin which T(x) is the transmission which depends on the intensity of the external influence: T(x) = xout/x. The appearance of such curves is largely determined by the total absorption coefficient k(x). Depending on the linear transmission T0, three types of curves can be observed corresponded to blacking-out (T0 < 0.5), limiting (0.5 < T0 < 0.8), and (non)linear (T0 > 0.8) attenuation of the input radiation. In the first and last cases, the absorber can be used as a radiation interrupter and a light filter, respectively.

Nonlinear Absorbance

Absorbance at a wavelength of 532 nm, being typical for providing laser experiments, demonstrates a nonlinear pattern concerning the strength of external electric fields (recall that we are considering the intensity of laser radiation, implying only its electrical component). The fitting of calculated absorbance with eq gives the smoothed curves (R2 > 0.9). Deviation from linearity is observed at high values of electric fields, which can be associated with an orbital distortion.[30,31] This, in turn, leads to disturbing the distribution of electron clouds involved in orbitals, which results in a higher overlap during electron transitions.[32] For this reason, in simulations, we have a picture similar to multiphoton absorption, which is observed in laser experiments.

Next, we need to cut off nonlinearity to follow the 2PA absorption model.[16,18] Differential analyses of eq have indicated two characteristic points on the absorption curves. The first, PNL, corresponds to the start of nonlinear processes, and the second corresponds to the saturation of absorption, PD:where t = ln (x).

The chart of the eq in linear coordinates is shown in Figure with a violet line.

Figure 2

Nonlinear absorption of dye 1c under different electric fields applied: TDDFT calculated spectral response is displayed with red dots; the blue and violet curves represent fittings corresponding to eqs and 7, respectively. The optical pathlength and linear transmission at 532 nm were set to 0.2 and 0.7 cm, respectively. The 2PA absorption model fit is shown with a green line (eq ). PD is the saturation of absorption, followed by bleaching, or optical breakdown of the absorber. Point PNL corresponds to the deviation of absorbance from linearity.

The extrema of d2k/dt2 are marked on the charts: maximum, PNL = exp (a – b), and minimum, PD = exp (a). Let us expand eq in a Taylor series in powers of n = 1 in the vicinity of the point PNL:

From the absorption curves, most parameters can be derived for assessing the effect of laser radiation on dyes (Table , eqs –14). However, there are still limiting thresholds that are also important for evaluating the effectiveness of optical limiters. We will search for them based on the analysis of the limiter curves.

Table 1

General Derived Parameters Based on the Applied Pure 2PA Absorption Model to Evaluate the Optical Limiting Properties of Dyes 1a–fa

parameter	equation (eq no.)
Calculated from absorption curves (k vs. input power)
linear absorption coefficient	9
nonlinear two-photon absorption coefficient	10
linear transmission (d is the optical pathlength)	11
optical destruction of the absorber	12
transmission at the PD	13
attenuation coefficient	14
Calculated from output curves (output vs. input power)
point of deviation from T0 on the output curve (limiting threshold)	15
point of completion of deviation from linearity	16
limiting threshold for the ideal optical limiter	17
dynamic range of the limiter	18
limiter activation speed	19

a, b, and c are the fitting parameters of eq .

Optical Limiting Effect

At high fields, dyes 1a–f show a drop in absorption observed a little before the saturation point PD (Figures and 3). We call this point the limiter destruction threshold. At P > PD, the limiter curves deviate upward, which in terms of a laser experiment represents the optical breakdown of the material. In Figure b,c, the area of the possible destruction of the material is shown with a dashed line. This result is partially consistent with the deviation from the applied pure 2PA absorption model (eq ) describing the overall process.

Figure 3

TDDFT simulation of some possible effects of the interaction of laser radiation with an absorber for different linear transmission values (an example for dye 1c): blacking-out (a), limiting (b), and nonlinear attenuation (c). The trajectory of an ideal limiter is shown with a brown line. Localization of points P0, P1, PD, and PL (Table ) was performed according to eqs and 15–17. The pure 2PA absorption model (eqs and 8) was utilized to fit output curves (blue solid lines). Deviation from the applied model is indicated with a broken line. The yellow background in (b) shows the area responsible for the speed of the limiter activation.

Destruction point PD (eq ) is important for calculating the dynamic range of optical limiters. However, first of all, it is necessary to find on the limiter curve the point of deviation (P0) of the transmission from the given linear value at the irradiation wavelength. To do this, we switch to logarithmic coordinates by the argument (x = e) and process a differential analysis. To simplify the expressions, we will work with eq . In logarithmic coordinates for the case of two-photon absorption, the limiter curve will be described by the following expression:

The derivatives of eq for t are represented as

That is, they are polynomials of the nth degree by the argument e( = x). Functions d/dt have extrema. The first of these is always the maximum, and it makes practical sense, while the others are not. The solution of equations for n > 2 is possible only numerically, and this gives the position of the potential point of deviation of the output signal from the linear transmission, i.e., limiting threshold (P0). The larger the n, the closer to zero this point tends. In general form, the solutions of these differential equations can be written as follows:where m is a real number (m ≤ 1). For m = 1, we have a position of point P1 (eq )—the completion of deviation from linearity, which is blurred in a certain interval of input energies. To localize the point P0 on the limiter curve, we introduce the parameter ξ, the absolute deviation of two functions: linear transmission (y = T0x) and nonlinear attenuation (eq ) at the point x = x0. The equationcannot be solved analytically for m, but the expression for the relative error δ looks much simpler:

Here, m → 0: δ → m. For m = 0.0760038 (the first extremum of the function d3y/dt3),

δ̅ ≈ 0.0760770. Such accuracy is quite enough for us to determine the deviation threshold P0, and this does not contradict visual control (see eq and Figure ). The ideal limiter curve is another indicator of the effectiveness of an absorber in optical limiting. Comparison with the real curve shape allows judging how quickly the limiter “turns on”. This is important for instantaneous protection operation. The limiting threshold for the ideal limiter, PL, is calculated by eq ; it is located at the intersection of two tangents marked in Figure 3 with brown lines. It is important to note that the thresholds P0 and PL are dependent on the nonlinear response of the material. The higher the β, the lower the thresholds, and therefore, at lower input radiation fluxes, the limiter starts to work. In turn, the rate at which the limiter is turned on (eq ) is directly proportional to β and also depends on linear absorption (i.e., concentration). However, a significant increase in α results in darkening (Figure a). Therefore, the issue of concentration and aggregation of dyes requires special attention along with their nonlinear properties.

The mathematical details regarding this section can be found in the Supporting Information. All analytical mathematical calculations were carried out with the Wolfram Mathematica 10.0 suite (https://wolfram.com/mathematica).

Results and Discussion

DFT Structures

The DFT optimized structures of dyes 1a–f are shown in Figure . As we noted in ref (17), slipped-cofacial phthalocyanine J-type dimers can exist as two isomers—parallel and oblique. These isomers belong to different symmetry groups; they are formed with equal probability in syntheses, but the composition of the mixtures cannot be easily controlled. The oblique isomer belongs to the C2 symmetry group, while the parallel isomer belongs to the C group. In the latter case, the structure has an inversion center, and therefore, second-order nonlinear optical properties are not realizable for it. Adjusting the parallel/oblique ratio in prospect will selectively control the magnitude of the second order nonlinear response. The isomerism of J-type cofacial dimers was first revealed by Kazuya et al. on the example of imidazolyl-substituted phthalocyanines.[33]

Structures of J-dimers 1a–f are characterized by the distortion of macrocycles, which prevents the aggregation of molecules in solutions and the solid phase. As a result, the ordered structures are formed.[16]

Optical Limiting Effect and Absorption Efficiency

In the present investigation, we assume that under the action of external fields, a molecule passes into excited metastable states, the lifetime of which is rather long compared to a light pulse. As the intensity of the external fields increases, the dipole-forbidden levels are populated in the visible region, and these electronic transitions contribute to the total absorption, enhancing it. It has been demonstrated that such an assumption leads to correctly reproducing the trend and the overall shape of the optical limiting fingerprint.[15] Placing a molecule in a strong electric field opens up exciting nonperturbative phenomena such as inconsistent multiple ionization, the origins of which lie in the subtle ways of interaction of electrons with each other. The potential of external electric fields is added to the Hamiltonian, and then the time-dependent Schrödinger equation is solved for the many-electron wavefunctions.[25] The applied electric fields induce charge transfer in the molecules, with the dipole moment of the molecule changing linearly at smaller fields and exponentially otherwise.[35]

The calculated TDDFT spectra of the model phthalocyanine J-dimers 1a–f under electric fields applied are shown in Figure . The response is presented in ΔA units: ΔA = A(F) – A(F = 0) to demonstrate the possible occurrence of the optical limiting effect (ΔA > 0), which is detected in a wide range of wavelengths. J-dimers 1a–f may be of interest as nonlinear absorbers in the 400–600 and 750–850 nm ranges where transparency is observed in the absence of incident radiation. At a wavelength of 532 nm (the second harmonic of an Nd: YAG laser), the absorption density for most dyes is high, which explains the good optical limiting response that we observed earlier.[18] However, nickel (1e) and cobalt (1f) complexes demonstrate low efficiency due to poor absorbance density at 532 nm.

Figure 5

TDDFT simulation (B97/6–31 + G(d)) of nonlinear absorbance (a) and optical limiting output signal (b) based on the calculated electroabsorption spectra (c) of phthalocyanine J-dimers 1a–f. The optical pathlength and linear transmittance at the analytical wavelength (532 nm) were set to 0.2 cm and 0.7, respectively.

With an increase in the strength of the electric field, the absorption response increases nonlinearly. However, in this study, we are restricting with a pure two-photon absorption model, which can be applied only for weak fields (0.001–0.005 a.u.). A deviation from this model allows finding the switch-on point (P1, eq ) on the limiter curve, after which an almost stationary response of the absorber to the external impact follows. Saturation absorption is the maximum on the absorption curve; further (at even stronger fields), a drop in absorption is possible, and this means an optical breakdown of the material (Figure c).

A more detailed consideration of the influence of electric fields on the nonlinear properties of molecules in TDDFT calculations will require the use of complicated models, in particular, (2 + 3) PA, for which it is necessary to add a quadratic term in eq or, say, removing linearity (β = 0), leaving only nonlinear terms. In our further works, we will pay attention to evaluating the contribution to limiting multiphoton processes, which can be simulated by using strong fields in calculating electronic transitions.

Correlations

Applying the factor analysis of the data presented in Table has allowed us to solve a multiparametric problem and simulate the relationship between the optical limiting parameters calculated by quantum-chemical modeling of nonlinear absorption and the ones obtained experimentally. Experimental and simulation data were used as dependent and independent variables, respectively. Employing a multi-iterative procedure of combining regressors implemented in the Correlato program,[36] strict correlations with experiments are found for the development of a forecasting model. Thus, having at their disposal only the structure of a substance, a researcher can figure out (1) whether this compound is capable of giving the optical limiting effect and (2) what the nonlinear response will be in practice. In other words, within the framework of a given model, with some approximation, it is possible to obtain experimental data without performing a costly experiments.

Table 2

General Computed Data for the Structures of Dimeric Phthalocyanines 1a–f and Experimental Data of the Nonlinear Properties of these Dyesa

	dyes
parameter	1a	1b	1c	1d	1e	1f
Optimized structures (B97-D/6-31 + G(d))
angle of slippage (deg)	24.1	21.2	22.4	23.8	36.5	20.8
tilt angle (deg)	3.9	17.1	13.2	4.1	1.1	12.9
relative rotation of the macrocycles (deg)	46.2	47.3	50.5	49.3	40.8	48.2
distance between macrocycles (Å)	4.8	3.8	3.9	3.8	3.8	3.9
Experimental data
average size of aggregates (nm)[34] (for substitution in eq 2)	53	124	105	116	131	70
radius of the laser beam in focus (μm)[18] (for substitution in eq 2)	64	67	64	68	48	64
nonlinear two-photon absorption coefficient (β, cm GW–1)[18]	315	340	360	228	57	21
limiting threshold (P0, GW cm–2)[18]	0.03	0.05	0.03	0.03	0.65	0.60
dynamic range of the limiter (DR)[18]	1000	460	830	930	72	82
attenuation coefficient (kA)[18]	8.0	7.0	7.8	6.4	3.2	1.7
simulationb
nonlinear two-photon absorption coefficient (β, cm GW–1), eq 10	8.09	10.31	12.12	6.93	3.01	2.38
optical destruction (PD, GW cm–2), eq 12	2.56	1.75	1.62	2.99	15.38	19.70
attenuation coefficient (kA), eq 14	8.41	7.71	8.01	7.33	4.77	3.65
limiting threshold (P0, GW cm–2), eq 15	0.04	0.04	0.03	0.05	0.25	0.21
dynamic range of the limiter (DR), eq 18	52.5	54.6	51.6	54.6	70.1	78.3
limiter activation speed (R, cm2 W–1), eq 19	11.0	13.3	16.2	9.6	2.2	2.3

Nd: YAG laser source; linear transmittance, 0.7 at 532 nm analytical wavelength; solvent, THF; concentration, ca. 2.3 × 10-4 mol L-1; optical pathlength, 0.2 cm.

Based on equations from Table .

Figure shows the result of the multivariate regression analysis. Thus, for the nonlinear absorption coefficient β and the radiation attenuation coefficient kA, the correlations are linear (R2 > 0.98), while for the limiting threshold P0 and the dynamic range DR, we failed to obtain a linear response. Perhaps this is due to the accuracy of their estimation in the experiment. In particular, to calculate the DR, it is necessary to know the exact value of the radiation energy destroying the limiter material, and it can be determined only approximately since the medium in which the absorber is placed also exhibits nonlinear properties at high input energies. The speed of switching on the limiter, r, is included in all correlation expressions (Figure ). The exclusion of this variable from the analysis greatly worsens the sought relationships. The highest values of r were found for zinc J-dimer 1c (16.2 cm2 W–1). The worst result was obtained for the nickel (1e) and cobalt (1f) complexes (∼2 cm2 W–1). This is due to the low degree of π–π interactions of macrocycles in these dimeric structures compared to other complexes. Although dyes 1e,f are more resistant to laser radiation, it is impractical to use them as optical limiters.

Figure 6

Experiment vs. simulation correlations for the nonlinear absorption coefficient β (cm GW–1), limiting threshold P0 (GW cm–2), dynamic range DR, and attenuation coefficient kA.

Conclusions

Thus, in this study, we have implemented an integrated approach aimed at quantifying the efficiency of optical limiting based on TDDFT calculations. Having only an optimized structure, first, the researcher, can reveal whether the nonlinear absorber is capable of demonstrating the optical limiting effect. Further, having chosen a certain standard, one can evaluate the relative efficiency of the optical limiting. However, for creating the forecasting model to obtain experimental values (with a certain degree of reliability) without providing experiments, it is still necessary to conduct laser experiments for a small series of similar absorbers under given conditions and build correlations. Obviously, the wider the series, the more accurate the prediction will be. Thus, the correlations between theory and experiment will provide the numerical result for some derivatives.

Numerical processing of TDDFT spectra is greatly simplified if a pure two-photon absorption model is used, resulting in cutting off the unnecessary nonlinearity in absorption caused by strong electric fields, especially since this nonlinearity leads to bleaching in simulation (Figure b,c) or optical breakdown in practice. However, using our approach, the researcher can derive similar equations, for example, for the case of pure three-photon or even multiphoton absorption, since the Taylor series expansion of absorption nonlinear function (in this work eq ) can be done with any number of terms.

We have implemented an indirect method for simulating the optical limiting effect, which is based on the calculation of absorption spectra with external electric fields that simulate the laser irradiation. This approach, firstly described by Cocchi et al.,(15) has the advantage that the researcher does not need to characterize all the phenomena occurring at the interaction of intense light with the substance; it is general and is not limited to phthalocyanines since TDDFT spectra can be calculated for each structure in the presence of electric fields. Also, it does not matter how many electrons a molecule has or whether it has an open or closed shell structure, since the TDDFT and finite field DFT methods are equally applicable to all structures. Only the calculation time can serve as a limitation since it is necessary to solve a large number of quantum-chemical tasks.