Process intensification in photocatalytic decomposition of formic acid over a tio2 catalyst by forced periodic modulation of concentration, temperature, flowrate and light intensity

The effect of forced periodic modulation of several input parameters on the rate of photocatalytic decomposition of formic acid over a TiO2 thin film catalyst has been investigated in a continuously stirred tank reactor. The kinetic model was adopted based on the literature and it includes acid adsorption, desorption steps, the formation of photocatalytic active sites and decomposition of the adsorbed species over the active titania sites. A reactor model was developed that describes mass balances of reactive species. The analysis of the reactor was performed with a computer-aided nonlinear frequency response method. Initially, the effect of amplitude and frequency of four input parameters (flowrate, acid concentration, temperature and light intensity) were studied. All single inputs provided only a minor improvement, which did not exceed 4%. However, a modulation of two input parameters, inlet flowrate and the acid molar fraction, considerably improved the acid conversion from 80 to 96%. This is equivalent to a factor of two increase in residence time at steady-state operation at the same temperature and acid concentration.


Introduction
The direct conversion of the small molecules (CO 2 , CH 4 , formic acid, N 2 , etc.) using renewable energy may be realized using three possible approaches: electrocatalysis, photocatalysis and non-thermal plasma (NTP), although a limited amount of data still exist for the comparison of all three alternative pathways in terms of conversion and energy efficiency. These three methodologies are often considered as complementary rather than alternative, and plasma and photocatalysis could be integrated, for example, to produce hydrogen in a distributed way. They are all part of the development of the so-called renewable energy-driven (or solar-driven) chemistry and energy [1]. Non-thermal plasma is an ionized gas created in mild conditions by applying electric energy to a gas. The plasma route is preferable for future industrial exploitability as it could be easily scaled up and it provides higher production rates and feed flexibility [2]. It is a potentially crucial technology to develop future chemistry based on the use of renewable energy and that is suitable for distributed use [3]. Photocatalysts are often investigated as proxy to NTP-driven catalysis as well as for benchmarking and for assessment of plasma routes. Titania thin films have become one of the most studied catalysts due to them being efficient, largely available and stable in a wide range of operating conditions including UV irradiation [4][5][6][7][8]. Specific novel thin film titania-based composite catalysts with embedded ferromagnetic nanodomains were also developed to induce radiofrequency heating and fast temperature modulations directly at the catalyst active sites [9][10][11]. By embedding ferromagnetic nanodomains in the catalyst support and inducing radiofrequency heating, the composite catalysts allow very fast heating and cooling rates in structured chemical reactors [12]. This allows for the periodic modulation of temperature in chemical reactors that could lead to significant improvement in reactant conversion and selectivity, thus introducing a further element to improve the catalytic performance over stationary conditions [13]. The modulation of input parameters was shown to be an effective tool to increase the performance and conversion in many applications [14][15][16].
The forced periodic operations could be realised in the form of a single or multiple parameter modulation(s). Zuyev et al. [13] developed a transient material and heat transfer reactor model to maximize the time-averaged production rate at minimum operational costs. Reactor modelling allowed optimization of the control strategy to improve the reactor performance in comparison to the steady-state operation. The reactor dynamics were described by a nonlinear system of ordinary differential equations controlled by two inputs: the inlet concentration and the inlet temperature in a chemical reaction of the type "A → product" [13]. The optimisation goal was to maximize the conversion of A to the product over a specified period of time. The problem was analysed by using Pontryagin's maximum principle with Lagrange multipliers. A control design scheme for maximizing the performance was proposed and verified by numerical solution of the reactor equations. This is an important aspect of process intensification [17,18]. This paper utilizes a nonlinear frequency response (NFR) method for analysis of the reactor performance in a reaction of photocatalytic decomposition of formic acid over a titania thin film catalyst. The NFR method is an approximate and analytical method, which is applicable only for stable, weakly nonlinear systems. In such systems, all nonlinear terms in the dynamic model equations can be represented in a polynomial form or can be expanded in the Taylor series form [19,20]. The theoretical basis and the principle of the nonlinear frequency response method was explained in our previous publications [15,19,21].
More recently, the so-called computer-aided nonlinear frequency response (cNFR) method was introduced by developing a user-friendly software application for implementing the NFR method [27]. By combining the NFR method and multi-objective optimisation techniques, a new methodology for the optimisation of periodic operations was established [36], which is much faster than the classical methods based on numerical integrations.
In this study, in order to evaluate the theoretical concept, we started with the simplest reactor applied in chemical reaction engineering, the continuously stirred tank reactor (CSTR). We assembled a kinetic model for photocatalytic oxidation of formic acid over a titania catalyst valid in the 20-130 • C temperature range. The selection of the reactor parameters was based on the possibility to manufacture the actual reactor system at Processes 2021, 9,2046 3 of 25 a laboratory scale. The performance of titania photocatalysts in this model reaction is quite typical for other photocatalytic (and plasma-catalytic) reactions and, thus, should allow the performance of some generalization with respect to the main trends of forced periodic operation.
Considered here is the effect of four input parameters: the inlet flowrate, the inlet concentration of formic acid, the process temperature and the intensity of UV light. The main contributions of this paper can be summarized as follows. In Section 2, we formulate the kinetic model of formic acid decomposition over titania thin films. A reactor model for non-stationary CSTR is described in Section 3. In Section 4, we give the main facts needed to understand the application of the NFR method, which is in the focus of this paper. In Section 5, the application of the NFR method for evaluating forced periodic operations of a CSTR reactor is presented. Section 6 describes forced periodic modulation analysis for selected cases (four cases with single input parameters and six cases with double input parameters).

Kinetic Model of Formic Acid Decomposition over a Titania Photocatalyst
The mechanism of formic acid decomposition over metal catalysts was reported in many studies and Tang et al. made a detailed analysis of the key factors responsible for high activity over bulk metal catalysts [37]. Considerably less attention was devoted to oxide catalysts, including titania. Oxidative dehydrogenation requires weak acidic properties [38]. Extensive investigations on rutile TiO 2 have demonstrated dissociative adsorption of formic acid, with formate species (HCOO) adsorbing in two types of bidentate bridging configurations. However, work on anatase TiO 2 surfaces is more limited, and conclusive results regarding the geometries and energetics of formic acid on anatase have not been reached, to the best of our knowledge. Vittadini et al. studied formic acid adsorption on both clean and hydrated anatase [39]. Their results indicate that the adsorption is purely molecular on the clean surface, while dissociation is favoured on the hydrated surface. In their subsequent study, they observed that formic acid dissociates very easily and adsorbs in different configurations onto the anatase TiO 2 (001) surface [40].
A simplified mechanism of the photocatalytic decomposition of formic acid consists of three steps (Plan 1) [41][42][43][44][45]. The first (reversible) step includes the interaction of formic acid from the gas phase, A, with adsorption sites on the photocatalyst surface, [ ], leading to the formation of formate species, [A], on these sites. In parallel the UV irradiation of surface sites converts them to photocatalytic active sites, z * (step 2). The last step consists of the decomposition of surface formates on the photocatalytic active sites accompanied by the evolution of products, carbon dioxide and hydrogen. The products are weakly bound to the surface and, once formed, are rapidly desorbed, resulting in the recovery of the active site. As the rate of step 2 is much higher compared to that of step 3, these two steps are often combined in a single step with an apparent kinetic constant.
Plan 1. Simplified mechanism of photocatalytic decomposition of formic acid. The adsorption isotherm is described by Equation (4) where ∆H is the enthalpy of adsorption of formic acid, and parameters P 1 and P 2 were obtained by fitting the temperature dependence of the equilibrium constant. The adsorption constant at 20 • C is equal to 1.34 × 10 −6 m 3 mol −1 s −1 , and the adsorption enthalpy was Processes 2021, 9, 2046 4 of 25 reported to be in the range 1.09-1.68 eV (105-162 kJ mol −1 ) [40]. As a reference value for our study, the minimum value from this range, corresponding to a bidentate bridging configuration, was chosen. The values of kinetic parameters are listed in Table 1. Popova et al. [46], by using in situ FTIR spectroscopy, observed formic acid and two types of formates (bidentate and unsymmetrical, Figure 1) during adsorption of 2.5 vol.% of formic acid on anatase with a surface area of 100 m 2 g −1 . They concluded that a single desorption step was not able to describe the desorption curve, therefore the desorption kinetics were described with two terms (Equation (5)).
Processes 2021, 9, x FOR PEER REVIEW 4 of 25 adsorption constant at 20 °C is equal to 1.34×10 −6 m 3 mol −1 s −1 , and the adsorption enthalpy was reported to be in the range 1.09-1.68 eV (105-162 kJ mol −1 ) [40]. As a reference value for our study, the minimum value from this range, corresponding to a bidentate bridging configuration, was chosen. The values of kinetic parameters are listed in Table 1. Popova et al. [46], by using in situ FTIR spectroscopy, observed formic acid and two types of formates (bidentate and unsymmetrical, Figure 1) during adsorption of 2.5 vol.% of formic acid on anatase with a surface area of 100 m 2 g −1 . They concluded that a single desorption step was not able to describe the desorption curve, therefore the desorption kinetics were described with two terms (Equation (5)).
The desorption of formic acid from an anatase catalyst was studied by temperatureprogrammed desorption (TPD) [47]. They identified two groups of desorption peaks. The first group appeared at 127 °C and was identified as formic acid. The position of this first peak allowed them to estimate the activation energy of desorption via a single point approximation for a first order desorption [48,49]  The desorption of formic acid from an anatase catalyst was studied by temperatureprogrammed desorption (TPD) [47]. They identified two groups of desorption peaks. The first group appeared at 127 • C and was identified as formic acid. The position of this Processes 2021, 9, 2046 5 of 25 first peak allowed them to estimate the activation energy of desorption via a single point approximation for a first order desorption [48,49] where k 0 is the pre-exponential term, θ is the surface coverage, β = dT dt is the heating rate and E A is the activation energy of desorption. At the maximum of the TPD peak, dr dt = 0, this corresponds to the following equation for a first order desorption: where T max is the temperature of the TPD peak. Assuming that the adsorption species are fixed on the surface and the pre-exponential term is equal to k D1 =1.0·10 13 s −1 , the activation energy of desorption (Ea D1 ) was estimated and listed in Table 1. The parameters of the second peak were found to be k D2 =2.05 s −1 and Ea D2 = 20920 J mol −1 .
Reaction rate data were also reported for formic acid decomposition over vanadia/titania catalysts (V 2 O 5 /TiO 2 ) with a high amount of titania, which were used for the rectification of kinetic data. Ivanov et al. studied the desorption of 1 vol.%. formic acid adsorbed onto V 2 O 5 /TiO 2 catalysts (80 wt.% TiO 2 ) with a surface area of 26 m 2 g −1 in the 50-150 • C range [50]. They reported desorption constants of 0.00523 and 0.01145 s −1 at 125 and 132 • C, respectively. These authors concluded that low temperatures are favourable, mainly for the recombination of bidentate formates with protons and formation of formic acid followed by its desorption, while the second peak is related to desorption from the bidentate formate species. The total desorption rate is determined by the sum of these two rates (Equation (8)). A factor 2 difference in the absolute values compared with the data reported by Popova et al. [46] can be explained by the presence of 20 wt.% vanadia in their catalysts.
Uslu [42] reported the value of the equilibrium constant for adsorption-desorption of formic acid over basic adsorbent, Amberlite IRA-67, a weakly basic gel-type polyacrylic resin with a tertiary amine functional group, which, in spite of a quite different chemical composition, has rather similar adsorption properties to those of titania. An adsorption enthalpy of −43,760 kJ mol −1 was reported. They also reported a pseudo first order adsorption constant value of 6.1·10 −4 s −1 at an acid concentration of 69.03 g L −1 (1.5 mol L −1 ). The adsorption kinetics were adequately fitted by the pseudo first order and pseudo second order adsorption models, however, the best result between the experimental and modelling data is achieved when using Equation (5).
The described kinetic parameters can be used to give the adsorption and desorption equilibrium rate equations r ads,1 = k ads (1 − θ A )n tot C A (8) where n tot is the density of adsorption sites on titania, where ρ is the site density, m is the photocatalyst weight in the reactor, A BET is the specific surface area, N A is the Avogadro number and V is the reaction volume. As an upper limit of the number of adsorption sites, the maximum value of the surface density of titanol groups Ti-OH, equal to ρ = 5.0·10 18 m −2 or 5 sites per nm 2 [51], was taken in this study. Upon UV light illumination, many oxide materials, including titania, can be used to decompose formic acid due to the photocatalytic effect via the pathway of charge carrier transfer [52]. The photoactivity of anatase and rutile was examined and interpreted by Sclafani and Herrmann [53] with reference to the densities of surface-adsorbed species. Their study demonstrated that higher levels of radicals adsorbed on the anatase surface gives rise to significantly higher photoactivity than rutile. Liu et al. [54] showed that the mechanism of formic acid photocatalysis over commercial P25 TiO 2 (primary particle size of >25 nm, surface area of~50 m 2 g −1 , anatase-to-rutile ratio of around 3:1) is the same in the low (30-100 • C) and elevated (100-150 • C) temperature range. The rate of thermal decomposition at 125 • C (4.0·10 −6 mol m −2 s −1 ) was almost two orders of magnitude below the rate of photochemical decomposition, therefore, the former can be neglected in the kinetic analysis. It was reported by several studies that, if a photocatalysis is driven by the charge carrier transfer, its rate depends on light intensity in a sublinear way according to where N is the exponent in the range typically between 0.2 and 1, k is a constant, I is the light intensity, and I 0 is the reference light intensity [55,56]. The rate of this process is therefore a property of the photocatalyst and independent of the reaction studied. Liu et al. [54] studied the photocatalytic decomposition of formic acid at 125 • C by varying the light intensity in the 1-5 mW m −2 range and found out that the index N was 0.3 in the whole range studied. This justified that photocatalysis is driven by the charge carrier transfer mechanism. The photocatalysis also follows the Arrhenius mode when the temperature was below 125 • C with an apparent activation energy of 8480 J mol −1 . These observations allow us to present the density of active photocatalytic sites formed in step 2 (Equation (2)) as follows.
Then, the rate of photocatalytic reaction (r 2 ) is modelled as a simple first order reaction, dependent on the recombination rate (k 2 ) and the density of photocatalytic sites (n R ) and the surface coverage (θ A ) of the photocatalyst with the formate species as follows It should be mentioned that, sometimes, the recombination rate is modelled as a second order reaction with respect to the concentration of holes [57]. However, the concentration of charges is not an appropriate measure since the charge carriers cannot freely move inside the reaction medium, but only inside their photocatalyst particles. Therefore, we chose to model the recombination rate as a function of the density of charges in the photocatalyst rather than a concentration of charge carriers.
Within the semiconductor, the intensity of light follows the exponential law [58] where I 0 is photon flux, α is the reciprocal absorption length, α = 1.0 × 10 4 cm −1 for titania [59] and L is the thickness of the photocatalyst layer, which can be varied in a rather wide range from a few nm films to micron thick particles. The photon flux depends on the frequency of the incoming light. For example, for a reference light intensity of 1.6 mW cm −2 at 365 nm, I 0 = 1.0 × 10 16 cm −2 . As the intensity reduces with the layer thickness, the magnitude of constant k 2 in Equation (12) would also depend on the configuration of the photocatalyst layer, and it may vary in a rather wide range between 0.001 and 0.24 s −1 [54]. The upper limit corresponds to a rather thin photocatalyst film of 5 nm.

Non-Stationary CSTR Model
As the dynamic reactor model provides the framework for the NFR method, it is important that for the temporal system, its assumptions are justifiable and correct. This section presents a system of equations describing the CSTR reactor coupled with the kinetics for formic acid decomposition outlined in Section 2, which gives a nonlinear Processes 2021, 9, 2046 7 of 25 system in relation to its inputs. The reactor volume (V) was taken as equal to 200 cm 3 , and the photocatalyst parameters for the titania film are listed in Table 2.
There are four possible modulated inputs, the inlet formic acid concentration, the inlet volumetric flowrate, the reactor temperature and the light intensity. This gives six binary combinations for simultaneous modulation of these four inputs: (i) inlet acid concentration and inlet flowrate; (ii) temperature and light intensity; (iii) temperature and concentration; (iv) temperature and flowrate; (v) concentration and light intensity; (vi) flowrate and light intensity. As the gas density changes with temperature according to the ideal gas law, the variation of flowrate and concentrations under periodic temperature modulation must be considered. Further considerations for a stochastic increase in gaseous flow are also required and are discussed later. Other assumptions include the same temperature and concentration in the reactor and in the outlet flow.
The overall dynamic mass balance within the bulk reactor, excluding adsorbate, is shown, Equation (15). Formic acid, H 2 , CO 2 and N 2 will be denoted with subscripts A, H, C and N, respectively.
Formic acid is the only adsorbate and, therefore, the overall surface mass balance is described by its adsorption ( where subscript 1 designates the adsorption-desorption steps (Equations (9) and (10), respectively) and subscript 2 designates the photocatalytic reaction step, i.e., the second step (Equation (13)). Similarly to Equation (16), the individual mass balances for all species in the reactor coupled with the kinetic rate equations, are described by Equations (17)- (24). . . . . .
The concentration fluctuations of formic acid with temperature are accounted for in the accumulation term, Equation (23). Substituting Equations (20)- (24) into Equations (15) Processes 2021, 9, 2046 8 of 25 and (16) gives the individual ith species mass balance in molar fraction terms, the species balances are generalised in Equation (25). The last term in Equation (25) gives gaseous expansion due to temperature changes and results from the derivation of Equation (23). Furthermore, this term can also account for changes in the gas volume where R A = r ads,1 + r des,1 , R H = r 2 , R C = r 2 , R N = 0. Then, the molar balances for individual components are given by Equations (27)- (30).
As discussed above, the outlet flowrate changes due to a stoichiometric increase (one mole of acid gives two moles of products, Equation (3)) and temperature modulation (accounted for in the end term in Equations (25), (27)- (30). Flowrate is derived by summing Equations (27)- (30) with an addition of the useful concept, ∑ y i = 1, which, when rearranged, gives the outlet flow term, accounting for increased and fluctuating gas flowrate.
Equations (26)- (30) and (32) describe the reactor behaviour and provide a basis for the computer-aided derivation of nonlinear FRF. Furthermore, these equations will be used to calculate the steady-state operating points for which the NFR method operates around.

The NFR Method for Evaluating Forced Periodic Operations
The method is based on the following facts: 1. When a stable system is perturbed by a periodic (sinusoidal or co-sinusoidal) input change (e.g., x(t)), after a long enough time the system output (e.g., y(t)) reaches a periodic (quasi) steady state [60]. If the system is weakly nonlinear, as are the majority of the processes in the chemical industry, the periodic steady state of the output is a sum of the basic harmonic of the same frequency as the input change and an indefinite number of higher harmonics and a non-periodic (the so-called DC) component [20] x(t) = x s + A cos(ωt) t→∞ → y(t) = y s + y DC + B I cos(ωt + ϕ I ) + B I I cos(2ωt + ϕ I I ) + . . . (33) where y s is the steady-state (SS) value of the output corresponding to the SS input value x s . A direct consequence of Equation (33) is that the time-averaged value of the periodic steady state of the output y is different from the corresponding SS value and that the difference is equal to the DC component Processes 2021, 9, 2046 9 of 25 This difference is a measure of the potential improvement owing to periodic input modulation. If the output of interest is, for example, the product flowrate at the outlet from a chemical reactor, the improvement would be possible only if y DC > 0. On the other hand, if the output of interest is, for example, the outlet reactant molar flowrate, the condition for improvement would be y DC < 0.
2. A convenient approach for the mathematical analysis of weakly nonlinear systems is the so-called concept of higher order frequency response functions (FRFs) (Weiner and Spina, 1980). This concept is based on replacing a nonlinear model of a weakly nonlinear system by a set of linear frequency response functions (FRFs) of different orders.
where G y,x, x, . . . , x n times stands for the n-th order FRF relating to the input x and output y.These FRFs are directly related to different harmonics and the DC component of the output defined in Equation (34) [20]. 3. For evaluating the possible performance improvement owing to the periodic modulation of an input, only the non-periodic (DC) component, which defines the timeaveraged value of the output, is of interest. The latter can be expressed as the following infinite series [15,19,22].
Equation (36) can be approximated with its first, most significant term [15,19,22] Equation (36) forms the basis for the nonlinear FRF method for fast evaluation of periodic operations with single input modulations. The function G (2) y,x,x (ω, −ω) is the socalled asymmetrical second order (ASO) FRF relating to input x and output y. This function can be derived from the nonlinear reactor model using a procedure described in [14,21,22]. The sign of the ASO function gives an answer as to whether periodic modulation can overperform SS operation. The magnitude of the possible improvement is proportional for the square of the input amplitude and can be estimated by Equation (37). 4. In the case of the simultaneous modulation of two inputs (e.g., x and z), the DC component of the output y is expressed as a sum of the contributions of the DC components related to the single inputs (x and z) separately, and the contribution of the DC component originating from the cross effect of both inputs [15,30] y DC = y DC,x + y DC,z + y DC,xz (38) It was shown that the best effects can be obtained if both inputs are modulated with the same frequency [15,30]. For co-sinusoidal modulations of inputs x and z, with frequency ω, input amplitudes A x and A z , respectively, and phase difference ϕ between them, these are described by Equation (39).
The separate contributions of the two inputs to the DC component can be approximately evaluated from the corresponding ASO FRFs, in the same way as in Equation (37).
Processes 2021, 9,2046 10 of 25 The contribution of the cross effect can be approximated by the corresponding cross ASO FRF [15,30] y,x,z (ω, −ω) (41) By choosing the appropriate value of the phase difference, it is always possible to obtain either a positive (when the output is to be increased) or a negative cross term (if vice versa). The optimal phase difference is calculated by Equation (42) for a positive value of y DC,xz or by Equation (43) for a negative value [30].
At the optimum phase difference, the DC component of the cross term increases with the increase in both input amplitudes. Combining Equations (38)-(41), an approximate expression for the DC component is obtained The DC component can be calculated for any chosen set of forcing parameters (frequency, amplitudes and phase difference) [24,30,31,34]. In case the output needs to be increased while one or both ASO FRFs corresponding to the separate inputs have negative signs, the input amplitudes also need to be optimised [15,21,30]. The optimal values of the phase difference and input amplitudes are functions of the frequency of the input modulations. In this work, the NFR method is used to analyse the potential for improving the performance of a continuous photocatalytic reactor by forced periodic modulations.

Application of the NFR Method for Evaluating Potential Forced Periodic Operations of a Photocatalytic Reactor for Formic Acid Decomposition
In this work, the NFR method explained in Section 4 is used for evaluating all possible cases of forced periodic operations with one or two modulated inputs, for the photocatalytic reactor for formic acid decomposition. The starting point for this analysis is the dynamic model presented in Sections 2 and 3. Defined in this section are all possible periodic operations, as well as the objective functions used for their evaluation. The goal of the analysis is to find the best solutions that would result in the highest conversion of the formic acid, taking into account the reactor capacity, as well.

Modulated Inputs and Outputs of Interest and Performance Criteria
For the CSTR, there are four possible inputs that can be modulated: the acid concentration, the volumetric flowrate, the temperature and the light intensity. For applying the nonlinear frequency response analysis, it is most convenient to define all inputs and outputs as dimensionless variables, as relative deviations from their corresponding SS values. For the case with four inputs, it is best to define them in vector form.
In Equation (45), all inputs are defined in the dimensionless form, as relative deviations from their steady-state values around which they are to be perturbed, which is more convenient for dynamic analysis in the frequency domain Based on the reactor model (Section 3), six outputs can be defined: (i-iv) the molar fractions of four components (formic acid, nitrogen, hydrogen and carbon dioxide) in the reactor, (v) the formic acid coverage and (vi) the outlet flowrate. However, as the photocatalytic step (Equations (2) and (3)) is irreversible, there is no need to evaluate the molar fractions of the products as they have no influence on the dynamics of the system. Therefore, the dynamic model can be reduced to three coupled equations (Equations (26), (27) and (32)) that allow us to evaluate the acid conversion, the acid molar fraction and the outlet flowrate. Thus, the other three outputs can be omitted from the analysis.
Furthermore, as the outlet flowrate changes with time, it is necessary to calculate its time-averaged value in order to obtain a correct value of acid conversion [21,24,25]. The formic acid outlet molar flowrate is in fact the main output of interest in our analysis, and it is defined by Equation (47).
Again, we define the outputs in a vector form The dimensionless form of the outputs, defined as relative deviations from their steady-state values, will be usedm

Possible Periodic Operations with One or Two Modulated Inputs
With four possible modulated inputs, it is possible to design a number of different periodic operations. In order to find the best possible scenario, all possible options with one or two modulated inputs will be investigated.
For single input modulations, there are four possible cases for periodic operation:

Frequency Response Functions for Evaluating Periodic Operations
For evaluating the formic acid conversion corresponding to any of the ten cases defined in Section 4, it is necessary to derive the asymmetrical second order FRFs corresponding to the dimensionless outlet acid molar flowratem A,out for single input modulations, as well as the corresponding cross FRFs, for Cases 5 to 10.
Nevertheless, because the model equations are coupled, in order to derive the FRFs corresponding tom A,out , it is necessary to derive the FRFs corresponding to the other outputs as well. Moreover, in order to derive the asymmetrical second order FRFs, it is necessary to first derive the first order FRFs, as the derivation procedure is recursive. Altogether, for our system with four possible modulated inputs and four outputs, it is necessary to derive:

•
Sixteen first order FRFs relating each output to each input; • Sixteen asymmetrical second order FRFs relating each output to each input; • Twenty-four cross asymmetrical second order FRFs relating each output to each combination of two inputs defined in Case 5 to Case 10.
The procedure for derivation of the first and asymmetrical second order FRFs is standard, and it has been described in detail in our previous publications [19,30,33,61]. In this work, the needed FRFs were derived using the cNFR software application, reported in [28]. These FRFs are defined to relate dimensionless outputs and dimensionless inputs.
Although fifty-six FRFs were derived, only the ASO FRFs corresponding tom A,out are needed for further analysis. The analytical expressions of the FRFs are rather cumbersome, so they are not presented here. Instead, the MATLAB function files, produced by the cNFR software, are provided in the Supplementary Information.
It is important to provide a reminder that the values of the ASO FRFs, which are functions of frequency, also depend on the steady state around which the system is perturbed and, naturally, of the physical and kinetic parameters of the system.

Performance Indicators
The performance indicator that can best express the efficiency of the investigated periodic operations is the acid conversion. For the forced periodic modulation, the conversion needs to be evaluated based on the mean (time-averaged) values of the outlet and inlet molar flowrates of formic acid, corresponding to the periodic steady state [21,24,25].
where . m A,in is the mean value of the acid inlet molar flowrate, defined as However, for Case 5, with simultaneous modulation of the inlet flowrate and molar fraction, it depends also on the forcing parameters As the inlet molar flowrate practically defines the reactor capacity, it will be used as an additional performance criterion. On the other hand, the mean value of the acid outlet molar flowrate where . m A,out,s is its steady-state value andm A,out,DC is the DC component of the dimensionless outlet molar flowrate. Using the NFR method, this DC component can be approximately evaluated based on the derived FRFs, in the following way: • For co-sinusoidal modulation of one input (x) with frequency ω and amplitude A x (Cases 1-4) • For simultaneous co-sinusoidal modulations of two inputs (x and z) with frequency ω, amplitudes A x and A z and phase difference ϕ (Cases 5-10) 1,x,z (ω, −ω) + sin(ϕ)Im G (2) 1,x,z (ω, −ω) , x = 1, 2 or 3, z = x + 1, x + 2 or x + 3, z ≤ 4 (56) By using Equation (54) with either Equation (55) or (56), the mean value of the acid outlet molar flowrate can approximately be calculated for any steady-state point and any set of forcing parameters (frequency, amplitude(s) and, for Cases 5-10, phase difference).
By introducing the expressions for . m A,in and . m A,out into Equation (50), the acid conversion can be calculated approximately.
For Cases 1-4 and 6-10 For Case 5 In Equations (56) and (57), x A,SS is the value of the acid conversion corresponding to the steady-state point around which the system is perturbed periodically.
As explained in Section 3, the NFR method is applicable only for stable systems. Consequently, it is necessary to check the system stability before applying the method. Along with the function files defining the FRFs, the cNFR also generates a file containing the system characteristic function [28]. The characteristic equation is obtained in the form of a linear quadratic equation. It is given in Appendix A.
The MATLAB function files for the ASO FRFs and cross ASO FRFs relating the outlet molar flowrate of formic acid to all four inputs and their combinations, generated by the cNFR software application, are used for simulations of the FRFs and evaluation of the formic acid conversion and optimisation of the forcing parameters. These results will be shown in the next section.

NFR Analysis for Different Cases
The NFR analysis was performed for a CSTR loaded with a TiO 2 thin film. The reactor operates at a pressure of 1 bar. The reactor and photocatalyst parameters are listed in Table 2.

Selection of the Steady-State Points for Analysis
Periodic modulation of the input variables always occurs around a steady-state point. Therefore, it is necessary to define several steady-state points for subsequent analysis. In this process, the steady-state temperature (T s ) was taken in the middle of the interval of the kinetic model (340 K), while the light intensity (I s ) was taken based on experimental values often cited in the literature (4.4 mW cm −2 ). Then, the acid molar fraction and the volumetric flowrate were optimised using multi-objective optimisation with two objective functions: the acid conversion (OF1) and the acid inlet molar flowrate (OF2). The SS optimisation was performed based on the SS reactor model, which is derived from the dynamic reactor model (Equations (27)-(30)) by setting all time derivatives equal to zero.
The multi-objective optimisation (MOO) of the steady-state operation was performed in MATLAB 2019b by using a non-dominated sorting genetic algorithm II (NSGA-II) [36,62]. The upper limit of y A,in,s parameter was limited based on the vapour pressure to prevent the formation of liquid phase at 340 K [63]. 0.010 < y A,in,s < 0.3189 The results of the steady-state multi-objective optimisation are presented in Figure 2, in the form of a Pareto front. Five points from this Pareto front (SS1 to SS5) were chosen for analysis of potential intensification by forced periodic modulation of the input parameters and are also shown as stars of different colours. The optimised SS values of the inputs (inlet acid molar fraction and inlet volumetric flowrate), as well as the resulting SS values of the outputs (the acid outlet molar fraction, the surface coverage and the outlet flowrate), together with the corresponding SS values of the performance indicators (the acid conversion and inlet molar flowrate) are listed in Table 3. The roots of the characteristic equation (Equation (A1)) are also listed to confirm the reactor stability in these points. As the characteristic roots are negative, the system is stable and the NFR method is applicable for the analysis of the reactor performance under forced periodic modulation for all steady-state points.

Evaluation of Periodic Modulations with a Single Input
As was discussed above, the sign of the asymmetrical second order FRF allows us to conclude on the reactor performance under periodic modulation. The corresponding FRFs are related to the acid outlet molar flowrate. Therefore, the periodic modulation will be superior to the steady-state modulation only if the corresponding ASO FRF is negative. Figure 3 shows the ASO FRFs corresponding to the four inputs defined in Section 5 (Cases 1 to 4), calculated for point SS5. It can be seen that the negative sign of the ASO FRF is obtained only for the case of temperature modulation (in the whole frequency range). The ASO FRFs corresponding to the other three inputs are positive for all frequencies. Consequently, only periodic temperature modulation may improve the reactor performance, while periodic modulations of other inputs would deteriorate the SS performance. Similar results were also obtained for the other four SS points (SS1, SS2, SS3 and SS4).

Evaluation of Periodic Modulations with a Single Input
As was discussed above, the sign of the asymmetrical second order FRF allows us to conclude on the reactor performance under periodic modulation. The corresponding FRFs are related to the acid outlet molar flowrate. Therefore, the periodic modulation will be superior to the steady-state modulation only if the corresponding ASO FRF is negative. Figure 3 shows the ASO FRFs corresponding to the four inputs defined in Section 5 (Cases 1 to 4), calculated for point SS5. It can be seen that the negative sign of the ASO FRF is obtained only for the case of temperature modulation (in the whole frequency range). The ASO FRFs corresponding to the other three inputs are positive for all frequencies. Consequently, only periodic temperature modulation may improve the reactor performance, while periodic modulations of other inputs would deteriorate the SS performance. Similar results were also obtained for the other four SS points (SS1, SS2, SS3 and SS4). Based on these results, only the periodic temperature modulation will be further analysed. The highest amplitude of the temperature modulation was set at 20 K (which can be experimentally reproduced). This corresponds to 5.88% of its SS value. Using this amplitude, the approximate DC component values for the outlet molar flowrate and the acid conversion were calculated using Equations (55) and (56), respectively, for the five chosen steady-state points ( Table 4). The relative change of the acid conversion versus the SS values and the frequency range corresponding to these conditions are also listed in Table 4.  Figure 4 compares the acid conversion at the 20 K temperature modulation versus the steady-state conversion around point SS5 in the frequency range between 10 −6 and 100 rad s −1 . It can be seen in Figure 4, that temperature modulation has very little impact on the acid conversion in the entire frequency range, although it is somewhat higher at higher frequencies. Similar results were obtained for the other four steady-state points (SS1 to SS4, not shown). It can be concluded that periodic modulation of a single input has little potential to improve the reactor performance. Therefore, the next section will explore the potential of simultaneous periodic modulation of two input parameters (Cases 5 to 10). Based on these results, only the periodic temperature modulation will be further analysed. The highest amplitude of the temperature modulation was set at 20 K (which can be experimentally reproduced). This corresponds to 5.88% of its SS value. Using this amplitude, the approximate DC component values for the outlet molar flowrate and the acid conversion were calculated using Equations (55) and (56), respectively, for the five chosen steady-state points ( Table 4). The relative change of the acid conversion versus the SS values and the frequency range corresponding to these conditions are also listed in Table 4.  Figure 4 compares the acid conversion at the 20 K temperature modulation versus the steady-state conversion around point SS5 in the frequency range between 10 −6 and 100 rad s −1 . It can be seen in Figure 4, that temperature modulation has very little impact on the acid conversion in the entire frequency range, although it is somewhat higher at higher frequencies. Similar results were obtained for the other four steady-state points (SS1 to SS4, not shown). It can be concluded that periodic modulation of a single input has little potential to improve the reactor performance. Therefore, the next section will explore the potential of simultaneous periodic modulation of two input parameters (Cases 5 to 10).

Periodic Modulation of Two Input Parameters
The analysis of Cases 5-10 is performed for periodic input modulations around points SS1 to SS5 to maximize the acid conversion. In the modulation of two input parameters, it is important to optimise the phase shift between these inputs in order to maximize the synergetic effect. It is equally important to optimise the amplitudes for both inputs as single input modulations could not improve the acid conversion. The optimisation was performed in the same frequency range from 10 −6 to 10 2 rad s −1 .
For Cases 6-10, the mean values of the inlet acid flowrate are equal to their SS values. Therefore, the phase shift corresponding to the maximum acid conversion is equal to that for the minimum acid outlet flowrate. The latter is calculated by Equation (43). However, the optimal amplitudes need to be obtained by a numerical solution. However for Case 5, the time-averaged acid inlet flowrate is not equal to the respective SS value as it depends on other forcing parameters (Equation (53)). Therefore, the optimal values of all three forcing parameters (phase shift and two input amplitudes) are numerically calculated and they are a function of the forcing frequency. The numerical solution was obtained by using the MATLAB function fmins. The following upper limits of the input amplitudes were used: 70% for the dimensionless inlet volumetric flowrate, 100% for the dimensionless inlet molar fraction, 5.88% for the dimensionless reactor temperature (corresponding to 20 K) and 100% for the dimensionless light intensity.
It was found that there exists no combinations of forcing parameters that would improve the process for modulations around steady-state points SS1 to SS5 for Cases 8-10, all involving light intensity as one of the modulated inputs. For the remaining cases, the optimisation results are shown in Figure 5. The highest acid conversion for the periodic steady state and the maximal gain versus steady-state operation (given in parenthesis) are shown. The corresponding frequency range is also identified.

Periodic Modulation of Two Input Parameters
The analysis of Cases 5-10 is performed for periodic input modulations around points SS1 to SS5 to maximize the acid conversion. In the modulation of two input parameters, it is important to optimise the phase shift between these inputs in order to maximize the synergetic effect. It is equally important to optimise the amplitudes for both inputs as single input modulations could not improve the acid conversion. The optimisation was performed in the same frequency range from 10 −6 to 10 2 rad s −1 .
For Cases 6-10, the mean values of the inlet acid flowrate are equal to their SS values. Therefore, the phase shift corresponding to the maximum acid conversion is equal to that for the minimum acid outlet flowrate. The latter is calculated by Equation (43). However, the optimal amplitudes need to be obtained by a numerical solution. However for Case 5, the time-averaged acid inlet flowrate is not equal to the respective SS value as it depends on other forcing parameters (Equation (53)). Therefore, the optimal values of all three forcing parameters (phase shift and two input amplitudes) are numerically calculated and they are a function of the forcing frequency. The numerical solution was obtained by using the MATLAB function fmins. The following upper limits of the input amplitudes were used: 70% for the dimensionless inlet volumetric flowrate, 100% for the dimensionless inlet molar fraction, 5.88% for the dimensionless reactor temperature (corresponding to 20 K) and 100% for the dimensionless light intensity.
It was found that there exists no combinations of forcing parameters that would improve the process for modulations around steady-state points SS1 to SS5 for Cases 8-10, all involving light intensity as one of the modulated inputs. For the remaining cases, the optimisation results are shown in Figure 5. The highest acid conversion for the periodic steady state and the maximal gain versus steady-state operation (given in parenthesis) are shown. The corresponding frequency range is also identified. Processes 2021, 9, x FOR PEER REVIEW 18 of 25 It can be seen that Case 5 shows a good potential for improvement as the acid conversion increases from 20 to 35%. On the contrary, a minor improvement was observed for Cases 6 and 7 in the entire frequency range. For that reason, only Case 5 will be considered in further analysis. Figure 6 shows the highest acid conversion at the simultaneous modulation of the acid molar fraction and inlet volumetric flowrate around point SS5. The steady-state conversion is also shown for comparison. It can be seen that the conversion is constant in the entire frequency range. The corresponding values of optimised amplitudes and phase shift are shown in Figure 7. The optimised amplitudes correspond to their highest possible values (100% for the acid molar fraction and 70% for the volumetric flowrate), while the optimal phase shift is very close to π rad (180 degrees), meaning that the modulations of the two inputs should be out-of-phase to achieve the best performance. Similar results are obtained for points SS1 to SS4 (not shown). The mean inlet acid flowrate for Case 5 is different from its steady-state value (Equation (54)). In fact, the mean value at the periodic modulation is lower than the steady-state value at a phase shift of π rad. For a proper comparison, the mean flowrates are also shown in Figure 7, together with their relative changes as compared to the respective steady-state values. It can be seen that Case 5 shows a good potential for improvement as the acid conversion increases from 20 to 35%. On the contrary, a minor improvement was observed for Cases 6 and 7 in the entire frequency range. For that reason, only Case 5 will be considered in further analysis. Figure 6 shows the highest acid conversion at the simultaneous modulation of the acid molar fraction and inlet volumetric flowrate around point SS5. The steady-state conversion is also shown for comparison. It can be seen that the conversion is constant in the entire frequency range. The corresponding values of optimised amplitudes and phase shift are shown in Figure 7. The optimised amplitudes correspond to their highest possible values (100% for the acid molar fraction and 70% for the volumetric flowrate), while the optimal phase shift is very close to π rad (180 degrees), meaning that the modulations of the two inputs should be out-of-phase to achieve the best performance. Similar results are obtained for points SS1 to SS4 (not shown). The mean inlet acid flowrate for Case 5 is different from its steady-state value (Equation (54)). In fact, the mean value at the periodic modulation is lower than the steady-state value at a phase shift of π rad. For a proper comparison, the mean flowrates are also shown in Figure 7, together with their relative changes as compared to the respective steady-state values.  The results of the periodic modulation of the volumetric flowrate and acid molar fraction are best illustrated on the Pareto front shown in Figure 2. Points PO1, PO2, PO3, PO4 and PO5, shown with circles of different colours, correspond to simultaneous periodic modulations of the inlet flowrate and the acid molar fraction around steady-state points SS1, SS2, SS3, SS4 and SS5, respectively. It can be seen that the simultaneous modulation of the flowrate and the molar fraction result in a substantial increase in conversion, although at the expense of production rate. Therefore, a multi-objective NFR procedure for optimisation of the periodic processes [27] is applied to optimise the production rate. In this procedure, two objective functions were optimized simultaneously by the selection of the steady-state point and the respective forcing parameters.

Multi-Objective Optimisation for Case 5
The multi-objective optimisation was performed using two objective functions: (i) the formic acid yield and (ii) the mean inlet molar flowrate. In this study, six variables were optimised. Two of them, the SS flowrate and the acid molar fraction, describe the position of the SS point, which should be used to achieve the highest improvement during periodic modulation. The other four are the forcing parameters (amplitude of the inlet flowrate,  The results of the periodic modulation of the volumetric flowrate and acid molar fraction are best illustrated on the Pareto front shown in Figure 2. Points PO1, PO2, PO3, PO4 and PO5, shown with circles of different colours, correspond to simultaneous periodic modulations of the inlet flowrate and the acid molar fraction around steady-state points SS1, SS2, SS3, SS4 and SS5, respectively. It can be seen that the simultaneous modulation of the flowrate and the molar fraction result in a substantial increase in conversion, although at the expense of production rate. Therefore, a multi-objective NFR procedure for optimisation of the periodic processes [27] is applied to optimise the production rate. In this procedure, two objective functions were optimized simultaneously by the selection of the steady-state point and the respective forcing parameters.

Multi-Objective Optimisation for Case 5
The multi-objective optimisation was performed using two objective functions: (i) the formic acid yield and (ii) the mean inlet molar flowrate. In this study, six variables were optimised. Two of them, the SS flowrate and the acid molar fraction, describe the position of the SS point, which should be used to achieve the highest improvement during periodic modulation. The other four are the forcing parameters (amplitude of the inlet flowrate, The results of the periodic modulation of the volumetric flowrate and acid molar fraction are best illustrated on the Pareto front shown in Figure 2. Points PO1, PO2, PO3, PO4 and PO5, shown with circles of different colours, correspond to simultaneous periodic modulations of the inlet flowrate and the acid molar fraction around steady-state points SS1, SS2, SS3, SS4 and SS5, respectively. It can be seen that the simultaneous modulation of the flowrate and the molar fraction result in a substantial increase in conversion, although at the expense of production rate. Therefore, a multi-objective NFR procedure for optimisation of the periodic processes [27] is applied to optimise the production rate. In this procedure, two objective functions were optimized simultaneously by the selection of the steady-state point and the respective forcing parameters.

Multi-Objective Optimisation for Case 5
The multi-objective optimisation was performed using two objective functions: (i) the formic acid yield and (ii) the mean inlet molar flowrate. In this study, six variables were optimised. Two of them, the SS flowrate and the acid molar fraction, describe the position of the SS point, which should be used to achieve the highest improvement during periodic modulation. The other four are the forcing parameters (amplitude of the inlet flowrate, amplitude of acid molar fraction, phase difference between the inputs and the frequency). The lower and upper bounds for all optimisation parameters are listed in Table 5.  Figure 8 shows the Pareto front for the multi-objective optimisation for Case 5. The five points on the Pareto front MO-PO1 to MO-PO5 (circles of different colours in Figure 8) were chosen for further analysis. These points have the same inlet molar flowrate as points SS1 to SS5 on the steady-state operation Pareto front (stars of different colours in Figure 8), which is also shown for comparison. It can be seen that a much higher acid conversion can be achieved by the forced periodic modulation of the acid molar flowrate below 2.5·10 −7 mol s −1 . The improvement is still possible at higher flowrates; however, the magnitude of this effect decreases.
Processes 2021, 9, x FOR PEER REVIEW 20 of 25 amplitude of acid molar fraction, phase difference between the inputs and the frequency). The lower and upper bounds for all optimisation parameters are listed in Table 5.   Figure 8 shows the Pareto front for the multi-objective optimisation for Case 5. The five points on the Pareto front MO-PO1 to MO-PO5 (circles of different colours in Figure  8) were chosen for further analysis. These points have the same inlet molar flowrate as points SS1 to SS5 on the steady-state operation Pareto front (stars of different colours in Figure 8), which is also shown for comparison. It can be seen that a much higher acid conversion can be achieved by the forced periodic modulation of the acid molar flowrate below 2.5·10 −7 mol s −1 . The improvement is still possible at higher flowrates; however, the magnitude of this effect decreases. The values of the optimisation parameters corresponding to the five selected points are shown in Figure 9. It can be seen that the outcome from the multi-objective optimisation gives somewhat different SS points than the single factor optimisation of the steadystate operation. The optimal input amplitudes are very close to their upper boundaries, while the optimal phase shift is close to 2 for all five points (meaning that the two inputs should be modulated practically in-phase). Figure 9 shows the conversion in points MO-PO1 to MO-PO5 and the relative difference from the respective SS values for the same inlet flowrates (points SS1 to SS5). It can be seen that a significant increase in acid conversion is possible at the same production rate. The conversion increases from 80 to 96% for The values of the optimisation parameters corresponding to the five selected points are shown in Figure 9. It can be seen that the outcome from the multi-objective optimisation gives somewhat different SS points than the single factor optimisation of the steady-state operation. The optimal input amplitudes are very close to their upper boundaries, while the optimal phase shift is close to 2π for all five points (meaning that the two inputs should be modulated practically in-phase). Figure 9 shows the conversion in points MO-PO1 to MO-PO5 and the relative difference from the respective SS values for the same inlet flowrates (points SS1 to SS5). It can be seen that a significant increase in acid conversion is possible at the same production rate. The conversion increases from 80 to 96% for point MO-PO1, and from 75 to 88% for point MO-PO2. The most important part of this result was achieved at a high conversion of acid when any further increase in conversion required a substantial (by a factor of 2) increase in residence time at steady-state operation to achieve the same level of conversion. This clearly shows a large potential of the forced periodic operation in process intensification of this highly important reaction. point MO-PO1, and from 75 to 88% for point MO-PO2. The most important part of this result was achieved at a high conversion of acid when any further increase in conversion required a substantial (by a factor of 2) increase in residence time at steady-state operation to achieve the same level of conversion. This clearly shows a large potential of the forced periodic operation in process intensification of this highly important reaction.

Conclusions
The performance of a CSTR for photocatalytic formic acid decomposition over a titania photocatalyst has been evaluated by forced periodic modulation of volumetric flowrate, acid molar fraction, light intensity and reactor temperature. A non-stationary dynamic reactor model has been developed that includes the photocatalytic reaction on the catalyst surface coupled with convective mass transfer and volumetric expansion terms due to the catalytic reaction. The frequency response functions were obtained by a computationally-aided nonlinear frequency response (cNFR) method. This allowed us to analyse the reactor response to the four inputs, either individually modulated or simultaneously modulated in pairs. The positions of the steady-state points of interest were further optimised using a multi-objective genetic algorithm to find the Pareto front of input variables that gives the highest acid conversion at the same flowrates as those during the respective steady-state operation. The individual modulation of either input parameter resulted in little or negative effect onto the acid conversion. However, the simultaneous modulation of inlet flowrate and acid molar fraction improved the acid conversion from 80 to 96%. This is equivalent to a factor of two increase in residence time at steady-state operation at the same temperature and acid concentration. The optimal values of forcing parameters (input amplitudes, frequency and phase shift) were identified. The optimal amplitudes were very close to their upper boundaries imposed by the technical limitations of the related control equipment, while the optimal phase shift was close to the complete 15

Conclusions
The performance of a CSTR for photocatalytic formic acid decomposition over a titania photocatalyst has been evaluated by forced periodic modulation of volumetric flowrate, acid molar fraction, light intensity and reactor temperature. A non-stationary dynamic reactor model has been developed that includes the photocatalytic reaction on the catalyst surface coupled with convective mass transfer and volumetric expansion terms due to the catalytic reaction. The frequency response functions were obtained by a computationallyaided nonlinear frequency response (cNFR) method. This allowed us to analyse the reactor response to the four inputs, either individually modulated or simultaneously modulated in pairs. The positions of the steady-state points of interest were further optimised using a multi-objective genetic algorithm to find the Pareto front of input variables that gives the highest acid conversion at the same flowrates as those during the respective steadystate operation. The individual modulation of either input parameter resulted in little or negative effect onto the acid conversion. However, the simultaneous modulation of inlet flowrate and acid molar fraction improved the acid conversion from 80 to 96%. This is equivalent to a factor of two increase in residence time at steady-state operation at the same temperature and acid concentration. The optimal values of forcing parameters (input amplitudes, frequency and phase shift) were identified. The optimal amplitudes were very close to their upper boundaries imposed by the technical limitations of the related control equipment, while the optimal phase shift was close to the complete angle (2π), meaning that the two inputs should be modulated in-phase. These results open up the avenue for further design of process control protocol. Additionally, the results demonstrate the power of the cNFR method to identify the optical case without conducting extensive simulations of the reactor performance at all possible combinations of input parameters and will serve as the basis for subsequent experimental verification.

Supplementary Materials:
The following are available online at https://www.mdpi.com/article/10.339 0/pr9112046/s1. The key data presented in this study are openly available in [wrap.warwick.ac.uk/160156] as MATLAB files. Folder S1_ FRFs-Provides the FRFs and includes a description of the contents and how to use. Folder S2_ Function files-Provides the MATLAB function files used to produce manuscript results and includes a description of folder contents on how to use.
Author Contributions: T.E. contributed to Sections 1-3 and 5-7. L.A.Ž. applied the cNFR method and performed numerical optimization (Section 6). P.D. contributed to Sections 1 and 2. R.S.A. contributed to Sections 1 and 3. M.P. contributed to Sections 1 and 4-6. E.V.R. contributed to Sections 1-3. All authors have read and agreed to the published version of the manuscript.

Funding:
The authors greatly appreciate the funding support from the Russian Science Foundation (grant number 20-69-46041). Thomas Ellwood and Evgeny Rebrov also acknowledge support from the ERC Synergy Grant "Surface-Confined fast modulated plasma for process and energy intensification" (SCOPE) from the European Commission with the Grant No. 810182.

Conflicts of Interest:
The authors declare no conflict of interest.