A Multiphysics Ray Optics Model for the Propagation of Electromagnetic Waves in Plasmas and the Design of Laser-Based Diagnostics in Nuclear Fusion Reactors

: Laser-based methods are widely used techniques for thermonuclear plasma diagnostics, since they can probe the internal of the plasma, being contactless and non-invasive. The interferometer, the polarimeter and Thomson scattering are the most widespread techniques, providing line-integral information of the electron density and the magnetic ﬁeld (interferometer–polarimeter) and local information of the electron density and temperature (Thomson scattering). The design of the diagnostics is a fundamental step, which usually requires an iterative process to maximise the performances of the diagnostics and their durability. In the future reactors, such as ITER and DEMO, the working environment will be much challenging, due to the various electro-mechanical, thermal and nuclear loads which may affect the life of the components and degrade the performances of the diagnostics. This work aims to present the modelling of plasma interferometry, polarimetry and Thomson scattering applied to a ray optics code. The model, developed on the COMSOL Multiphysics software, can be easily interfaced with multiphysics problems, allowing the possibility to test the performances of the diagnostics in several conditions.


Introduction
Plasma diagnostics is one of the most important research fields in nuclear fusion, since they are the inputs of the control system, which has the scope to control the confinement stability of the plasma, and they are fundamental to understand and analyse the plasma physics and equilibrium [1][2][3][4][5]. Contact measurements are prohibited for core plasma investigations, because they are invasive, their life would be too short, and their interactions with plasma would produce impurities, causing off-normal events, such as disruptions. Magnetic diagnostics provide only boundary or integrated information, causing a lack of knowledge of the interior of the plasma, making optical diagnostics one of the best methods to probe the plasma interior [6,7]. For example, bolometers are widely used to measure the total emission of the plasma [8]; fundamental for the power balance, and, by tomographic inversion methods, they are able to evaluate the radiation pattern on the poloidal plane [9][10][11]. Between the optical diagnostics, the laser-based ones found many applications, where the most valuable are interferometry, polarimetry and Thomson scattering (TS) diagnostics [12].
The interferometer, measuring the phase difference between the probing beam and the reference beam, is used to measure the line-integrated electron density along the line of sight. It is considered to be one of the most important plasma diagnostics [13,14]. For example, the Joint European Torus (JET) mounts four vertical and four horizontal interferometric lines of sight, which are used to probe the electron density of the plasma and control the plasma discharge. Plasma polarimetry is based on the birefringent characteristic of magnetised plasmas, which involves two components of the electric field of the probing beam propagate differently, involving a changing of the polarisation of the beam. Measuring the polarisation variation effects (Faraday rotation and Cotton-Mouton phase shift) of the beam, it is possible to obtain information about the line integrated electron density and magnetic field [15][16][17]. Moreover, the polarimeter and the interferometer can be performed on the same probing beam, as done at JET, making the interferometer-polarimeter a unique, reliable, and informative diagnostic. At ITER, it is expected to have different polarimeter-interferometer diagnostics, both toroidal and poloidal, and both for the study and the control of the plasma. The Thomson scattering analyses the scattering radiation of a probing beam. The scattered radiation intensity is directly proportional to the local electron density, while the scattered spectrum is dependent on the local electron temperature, because of the Doppler effect which occurs due to the interaction of the electromagnetic radiation with the moving electrons [18,19].
The design of these diagnostics is a crucial point of the building process of a tokamak since several factors must be considered. At first, since the interaction between the probing radiation and the plasma are strongly influenced by the plasma properties, it is important to simulate the effects at different plasma configuration and find the best setup which is optimised for all scenarios. For example, the laser wavelength of the polarimeter at JET is 195 µm, which is unusable for ITER or DEMO plasma scenarios [20][21][22]. Moreover, in the present tokamaks, the environment is not critical in terms of neutron and heat fluxes, as well as electromagnetic forces, while in the future tokamaks, these loads will be much higher and they will play an important role on the design of the diagnostics, since they may involve the degradation of the performances. The high neutron flux may change the optical properties of the mirrors, windows, and other optics, while the heat fluxes and the electromagnetic forces may vary the alignment, involving a decrease of the performances of the diagnostics [23][24][25]. Thus, a preventive effort to find more reliable and resilient solutions is fundamental to ensure the longest duration of the best performances. For these purposes, numerical simulations are decisive in finding the best diagnostic design. Given the various physics which need to be considered during the design of the diagnostics (optics, thermodynamics, neutrons, electromechanics, etc.), a numerical model able to easily implement multiphysics simulations will help in the development of the diagnostics of the future tokamaks.
This work aims to present the methodology implemented on the ray optics module of COMSOL Multiphysics to simulate three of the most important laser-based diagnostics in nuclear fusion reactors: interferometry, polarimetry, and Thomson scattering. The solutions presented and employed in COMSOL are general, and the same methods can be easily implemented in other ray optics models. The next section introduces the logic of COMSOL, and the equations used to simulate the various effects, while Section 3 analyses some numerical simulations with the aim to validate the logic implemented in the model. The validation will be performed by comparing the expected (or imposed value) with the measured one by the numerical simulation. At last, the most important issues and results are resumed in the last section.

Modelling
This section discusses the numerical models used to simulate the propagation of electromagnetic waves in a plasma. The numerical model has been implemented in COMSOL Multiphysics, a numerical simulation software which allows an easy interface of different physics. The model has been built over a standard COMSOL physics, the ray optics, and it has been implemented to take into account the propagation of electromagnetic properties in a thermonuclear plasma.

Ray Optics
The propagation of electromagnetic waves in a medium can be approached in two different ways, by the ray optics and the wave optics.
In the case of wave optics, or physical optics, the light propagates as a wave. Wave optics allows simulating, without modelling approximations, the exact propagation of electric and magnetic fields of light, and the major part of the optical phenomena, such as interference and diffraction, which normally cannot be simulated by geometrical optics. On the other hand, numerical models based on the wave optics are computationally demanding, since the maximum mesh element size must be comparable with the light wavelength. This issue limits the wave optics numerical simulations to applications with geometries that have sizes comparable with the light wavelength, making this approach unsuitable for thermonuclear diagnostics, where the ratio between the geometry size and the radiation wavelength (L/λ) is of the order of 10 7 for the TS, and a volumetric ratio (L/λ) 3 of 10 21 (considering only the vacuum vessel) [26].
Ray optics, also known as geometrical optics, describes the propagation of light as rays. The propagation of these rays is based on the reflection, refraction, and absorption equations. Ray optics is the common approach used to simulate the propagation of electromagnetic waves when the geometry size is orders of magnitude larger than the light wavelength, since the computational load required by the numerical simulation is much lower. However, without proper modelling, some physical effects, such as interference, diffraction and polarisation, cannot be simulated by ray optics [27].
The numerical model developed and presented in this work is based on a ray optics approach, where any effect due to wave propagation and interaction with the plasma is simulated as described in the following sections. The base model of ray optics simulates the propagation of the electromagnetic wave by solving the equations of the position (q) and wave vector (k) for each ray [27]: where ω is the angular frequency of the light and the local speed of rays is defined by the speed of light (c) and the local refractive index of the medium (n), i.e., v = c/n. Thus, the propagation of the ray position and wave vector is directly simulated once the refractive index field and the light wavelength (or frequency) are known. A magnetic plasma is birefringent and the refractive indexes of the "slow" (n 1 ) and "fast" (n 2 ) waves can be described by the Appleton-Hartree Equation (the collision term is neglected) [6,15,28] as follows: where ω p is the plasma frequency and ω c is the electron cyclotron frequency: And N e is the electron density, e and m e are the charge and mass of the electron, ε 0 is the dielectric permittivity and B 0 is the module of the magnetic field. F is calculated as: And θ is the angle between the magnetic field and the light propagation direction. In thermonuclear plasmas, the effect of magnetic field on the ray propagation can be neglected (to a first approximation) and the refractive index becomes [6]: The ray optics model is based on the continuous interaction between the rays and the "environment", i.e., the plasma, the vacuum, the air, the optics, etc. The ray properties are stored in a vector which evolves according with the propagation of the ray itself, while the environment properties are stored in the mesh of the numerical simulation. At the beginning of the simulation, the initial status of the ray at the time t 0 is known and the successive instant is calculated numerically with a specific numerical scheme (such as the backward Euler method). The solution of the ray optics equations requires both properties belonging to the ray and to the environment (such as the refractive index). The ray properties, being stored in the ray vector, are known at each time step, while the environment properties are calculated through the position of the ray at each time step.
The numerical simulations in COMSOL are principally divided into parameters, variables, geometry, domains, boundaries, physics and materials. The parameters are the constants of the simulation, which do not change during the simulation while the variables are solved during the numerical simulations at each time step. The geometry is used to design the environment, and it can be divided into domains. Each domain may have unique characteristics, such as plasma, air, vacuum, etc. The boundaries are used to simulate discontinuities, such as a wall, mirror, material change, etc. The physics defines the model equations and different physics can be used together for multiphysics purposes. For example, the mechanical deformation of a mirror due to a heat load on the structure can be simulated, and the variation of the ray reflection can be simulated as well, using mechanical, thermal and geometrical optics physics. Materials are used as libraries of material properties which can be assigned to the domains.

Interferometry
The interferometer is based on measuring the refractive index of the plasma by probing it through a laser beam. The phase of the electric field of the beam is [6]: and by comparing it with the phase of a reference beam, the phase shift can be directly linked with the plasma electron density [6]: where N c is called the cut-off density, defined as N c = ω 2 m e ε 0 /e 2 . When N e is much smaller than N c , and after taking into account the series development of the binomial, the equation can be approximated as [6]: If the electron density is larger than the cut-off density, the probing beam is strongly absorbed by the plasma and the interferometer ceases to work (n 2 is negative, thus n is imaginary).
In a wave optics approach, the base equations of the model would automatically compute the phase shift, while in the ray optics model, the phase shift must be modelled. In our model, an auxiliary variable, called φ, which computes the phase of each ray has been introduced. This variable is computed using Equation (8), where the refractive index of the plasma has been written as Equation (9). Thus, for any time instant t i , the position of the ray is x', y' and z' in the laboratory frame of the reference are known and the refractive index of the propagation medium (plasma, air, vacuum, glass, etc.) is calculated. The integration of Equation (8) is performed considering that dz is vdt, where dt is the time step.

Polarimetry
While the presence of the magnetic field can be neglected to evaluate the propagation of the ray, it is important for the evolution of the polarisation of the electromagnetic wave. In fact, the presence of the two refractive indexes, even if similar, changes the polarisation of the beam crossing the plasma. Plasma polarimetry is based on measuring the change of the polarisation to obtain some integral parameters directly correlated with the electron density and the magnetic field.
The Stokes notation has been used to model the polarisation variation in plasma [15,29]. The reduced Stokes vector has three components that are directly correlated to four polarisation angles: the phase shift φ, the electric field component ratio tan(α), the ellipticity tan(χ) and the ellipse orientation ψ, as follows [15]: The evolution of the Stokes vector in a magnetised plasma can be evaluated by the following differential vector equation [15]: The vector Ω takes into account the properties of the plasma and magnetic field [15]: where B p is the component of the magnetic field along the beam propagation direction, while B x and B y are the orthogonal components, the constant C 1 is 2.45 × 10 −11 and C 3 is 5.26 × 10 −13 [15]. The Faraday rotation (∆ψ) and the Cotton-Mouton phase shift (∆φ) can be directly correlated to the line integral of the electron density and magnetic field. The most common approximation, known as type 1, links the polarisation effects with the plasma parameters as follows [17]: The propagation of the Stokes vector in the model has been considered by means of three new auxiliary variables, one for each reduced Stokes component (s 1 , s 2 and s 3 ), using Equation (12). Attention must be paid to the reference frame used to describe the Stokes vector and the magnetic field. In nuclear fusion plasmas, the magnetic field is usually given in terms of a toroidal field (B t ), radial field (B r ) and vertical field (B z ), where r is the horizontal direction and z is the vertical direction on the poloidal section, while the Stokes propagation, described by Equation (12), is written in the reference frame of the probing beam. The direction of propagation of each beam is described by its unit vector as: One of the two orthogonal unit vectors, defined as the "x" direction, is calculated as n x =   n p,t −n p,r n p,z   , which is always orthogonal to the n p unit vector. The last vector, which describes the "y" direction, is calculated by the cross-product between n p and n x . Thus, the magnetic field is rotated from the reference frame of the reactor to the reference frame of the electromagnetic wave. Given the magnetic field on the reference frame of the electromagnetic wave, the electron density field and the laser wavelength, the Ω is calculated. The reduced Stokes auxiliary variables are calculated by an explicit backward Euler numerical scheme:

Thomson Scattering (TS)
Thomson scattering consists in measuring the scattered radiation of a laser beam. The electromagnetic wave, interacting with the charged particles, forces them to oscillate accordingly to the electromagnetic field, causing the emission of radiation at the same frequency of the incident wave. The thermal motion of an electron in a thermonuclear plasma makes the scattered spectrum broader than the incident radiation spectrum due to the Doppler effect. Moreover, due to the relativistic aberration (or headlight effect), the spectrum is blue shifted.
The spectral power, P s (R, ω s ), of the Thomson scattering radiation due to electrons can be written as follows [19]: where ω s is the frequency of the scattered radiation, P i is the power of the incident radiation, r 0 is the classical radius of the electron, L is a typical length associated with the scattering volume, N e is the electron density, dΩ is the scattering solid angle, ω is the difference between the scattered and the incident radiation frequencies, s is the unit vector of the scattered radiation direction, E i0 is the electric field vector of the incident wavelength, and k is equal to 2k i sin(θ/2), where k i is the wavenumber of the incident wavelength and θ is the scattering angle. The variable a is the mode velocity of the electrons in the scattering volume, defined as a = (2k b T e /m e ) 1/2 , where k b is the Boltzmann constant, T e is the electron density and m e is the electron rest mass. In practical application, the scattering spectrum is usually written in terms of radiation wavelength, using the scattered wavelength λ s and incident wavelength λ i , with ∆λ = λ s − λ i . Equation (18), for probing a beam linearly polarised at ϕ = 0, becomes [19]: Appl. Sci. 2021, 11, 434 The corr variable is used to consider the effect of the relativistic aberration. The complete correction, describe by Matoba et al. [29], can be written as: The numerical simulation of Thomson scattering is based on performing two steps. The first step simulates the propagation of the incident radiation while the second step simulates the scattered one. In the first step, four different accumulated variables are computed. These variables are the output of the first step and the input of the second step, allowing to determine all the characteristics of the scattered radiation. From Equation (19), it can be observed that the scattered spectrum is a function of the incident power, wavelength and scattering angle. The incident wavelength is defined as an input of the model. Three accumulated variables measure respectively the t, r and z directions (reference frame of the tokamak). The fourth accumulated variable measures the number of rays which crosses the mesh element (and thus the amount of incident power).
The second simulation releases the rays from the scattering volume with properties which are defined by the local properties of the plasma, the accumulated variables of the first simulation and the direction of the released rays. In order to have a spectral measurement, it is needed to use several rays at different wavelengths. Given a chosen wavelength, and the elements of the solid angle, a ray has the average power contained in these bins: where the index i defines the i-th direction and the index j the j-th wavelength element. A small number of the elements does not involve errors that are very large. Consider that in tokamaks the spectrum measurements are performed by a limited amount of spectral channels (from 3 to 5) and the temperature fit is obtained with a very small spectral resolution [18]. In practical applications, the electron temperature is calculated by fitting the spectrum data with the theoretical curves at different temperature and the best fit selects the temperature measured. Moreover, by measuring the total scattered intensities, it is possible to measure the local electron density, being known the scattering cross-section of one electron. However, because of the effect of optics, which may have a propagation spectrum that is not perfectly wide, the TS effect is usually calibrated [30].
The scattered radiation also has a change of polarization with respect to the input polarisation. This effect is strictly proportional to plasma temperature and a polarimetric TS could be used to measure the electron temperature of the plasma. However, the depolarisation of the TS radiation is a relativistic effect and it has not yet been implemented, since the temperature of the plasma obtained in present Tokamaks (relativistic depolarisation becomes interesting for electron temperatures higher than 15 keV.) The polarisation of the scattered radiation has been modelled considering the relativistic effects described in the literature [31][32][33][34]. Given the incident Stokes vector (S in ) and the plasma TS Mueller matrix (M TS ), the scattered Stokes vector is [34]: where the TS Mueller matrix is defined as [34]: where u = cos θ, µ = 1 τ = m e c 2 /T e and G(µ) = K 1 (µ)/(µK 2 (µ)), where K 1 and K 2 are the modified Bessel functions of the second kind. The definition of the scattered Stokes vector needs the accumulation of three other variables from the first step: the Stokes vector components (S in,0 is the total energy, yet calculated from the ray density). Please note that this is the frequency integrated Mueller matrix, which does not take into account the different polarisations at different wavelengths. This problem could be solved just using the frequency-resolved Mueller matrix [33]. However, theoretical studies showed that, for wavelength relevant for nuclear fusion, the polarisation difference is very small and completely negligible for the measurement purposes. Once the incident Stokes vector and the local properties of the plasma are known, the scattered Stokes vector is calculated accordingly with Equation (22). The propagation of the polarisation vector is computed through the model in Section 2.3.

Interactions with the Boundaries and Other Elements
The interaction of the ray with boundaries is governed by specific equations which depend on the type of boundaries. The most common and relevant for this work are walls, material discontinuities, mirrors, inlets and ray detectors. Walls have different possibilities. The ray may freeze, stick, disappear, have a specular reflection, have a diffuse scattering, pass through, etc. The disappear function is used to remove from the simulation the rays which touch the walls. The freeze function retains the ray position and direction (since they no longer change over time), while stick retains the ray position but resets all other variables to zero. The other wall functions are general cases of the following more specific boundary conditions. Material discontinuities are used to simulate Snell's law when a ray passes from one material to another. In this case, it is possible to require that both the refracted and reflected rays are generated (the incident ray becomes the refracted while the reflected is a new ray, called secondary). Mirrors specularly reflect the rays at the boundary. Inlets allow releasing rays with different properties. It is possible to set any ray feature, such as release time, the number of rays, initial direction, ray frequency and wavelength, initial phase, initial intensity and polarisation, and the initial values of the auxiliary variables. The ray detector returns a Boolean variable which describes if the ray has reached the ray detector, mostly used for post-processing. The boundary condition which involves a continuation of the ray propagation, such as material discontinuities and mirrors, allows defining the value of auxiliary variables after the interaction with the boundary condition. For example, the reduced Stokes vector components after the interaction with a mirror are defined by the incident Stokes vector components and the Mueller matrix of the mirror.

Validation on a Case Study
The models described in the previous section are used here in a simple case, to demonstrate the functioning and validate the model by comparing the model results with the expected ones. Three different numerical simulations have been performed, one for diagnostics.

Geometry, Domain, Boundary Conditions and Mesh
The geometry of the numerical simulations consists of an ITER-like vessel. The major radius of the Tokamak is 6.2 m, while the minor radius is 2 m [35]. The internal vessel has been built according to the poloidal section of ITER, while the external structure has been placed just to represent a solid approach and it does not play any role in the numerical simulation. The geometry of the numerical simulation is shown in Figure 1.

Validation on a Case Study
The models described in the previous section are used here in a simple case, to demonstrate the functioning and validate the model by comparing the model results with the expected ones. Three different numerical simulations have been performed, one for diagnostics.

Geometry, Domain, Boundary Conditions and Mesh
The geometry of the numerical simulations consists of an ITER-like vessel. The major radius of the Tokamak is 6.2 m, while the minor radius is 2 m [35]. The internal vessel has been built according to the poloidal section of ITER, while the external structure has been placed just to represent a solid approach and it does not play any role in the numerical simulation. The geometry of the numerical simulation is shown in Figure 1. The plasma used for the numerical simulation has been built applying a circular equilibrium type, widely used in the literature for the investigation of polarimetry performances. The electron density and temperature of the plasma have a parabolic profile. Defined by the dimensionless minor radius ( , ) = ( + ) . / , where r is the minor radius coordinate and a is the minor radius of the plasma, the electron density can be written as ( , ) = (1 − ( , ) ) and the electron density is ( , ) = (1 − ( , ) ) , where N0 is the electron density on the magnetic axis of the plasma. The toroidal field, defined by its value B0 on the major radius R0, has the form ( ) = . The profile of the plasma current density has a parabolic profile as well and it is calculated from the plasma current (Ip) as ( , ) = (1 − ( , ) ). The poloidal field is calculated from the integration of the current density and it is ( , ) = ( , ) ( , ) . Then, the refractive index field is computed (see Section 2, Equation (7)). Figure 2 shows the field of the electron density, toroidal and poloidal magnetic field, and of the differential refractive index, defined as the vacuum refractive index minus the plasma refractive index. The plasma used for the numerical simulation has been built applying a circular equilibrium type, widely used in the literature for the investigation of polarimetry performances. The electron density and temperature of the plasma have a parabolic profile. Defined by the dimensionless minor radius ρ(r, z) = r 2 + z 2 0.5 /a, where r is the minor radius coordinate and a is the minor radius of the plasma, the electron density can be written as N e (r, z) = N 0 1 − ρ(r, z) 2 and the electron density is T e (r, z) = T 0 1 − ρ(r, z) 2 , where N 0 is the electron density on the magnetic axis of the plasma. The toroidal field, defined by its value B 0 on the major radius R 0 , has the form B t (R) = B 0 R 0 R . The profile of the plasma current density has a parabolic profile as well and it is calculated from the plasma current (I p ) as j(r, z) = 2I p πa 2 1 − ρ(r, z) 2 . The poloidal field is calculated from the integration of the current density and it is B p (r, z) = µ 0 I p 2πa . Then, the refractive index field is computed (see Section 2, Equation (7)). Figure 2 shows the field of the electron density, toroidal and poloidal magnetic field, and of the differential refractive index, defined as the vacuum refractive index minus the plasma refractive index.

Interferometer Simulations
The numerical simulations of the interferometer have been performed by pla three vertical lines of sight, as shown in Figure 3. Three lasers emit beams with a w length of 20 µm and each laser emits 100 rays. Each beam interacts with a 50-50 b splitter. Thus, half of the rays are reflected, while the others are transmitted. The tran

Interferometer Simulations
The numerical simulations of the interferometer have been performed by placing three vertical lines of sight, as shown in Figure 3. Three lasers emit beams with a wavelength of 20 µm and each laser emits 100 rays. Each beam interacts with a 50-50 beam-splitter. Thus, half of the rays are reflected, while the others are transmitted. The transmitted rays interact with the mirrors on the top of the tokamak and enter in the vessel, probe the plasma, exit from the vessel and are reflected by the mirrors placed below the vessel. Then, they pass through the beam mixers and are collected by the detectors. The reflected rays from the beam-splitters represent the "reference" beams, used to measure the differential phase. They are reflected by the mixers and collected by the detectors.

Interferometer Simulations
The numerical simulations of the interferometer have been performed by placing three vertical lines of sight, as shown in Figure 3. Three lasers emit beams with a wavelength of 20 µm and each laser emits 100 rays. Each beam interacts with a 50-50 beamsplitter. Thus, half of the rays are reflected, while the others are transmitted. The transmitted rays interact with the mirrors on the top of the tokamak and enter in the vessel, probe the plasma, exit from the vessel and are reflected by the mirrors placed below the vessel. Then, they pass through the beam mixers and are collected by the detectors. The reflected rays from the beam-splitters represent the "reference" beams, used to measure the differential phase. They are reflected by the mixers and collected by the detectors. Nine numerical simulations were performed changing the electron density to simulate different plasma conditions, from the absence of the plasma to the flat top condition. The electron density (N0) goes from 0 to 2×10 20 m −3 . Figure 4 shows the results of the numerical simulations. The figure on the left represents the ray trajectories from the laser to the detector, where the ray colour shows the calculated phase according to Equation (8). Since the "reference" and the "plasma" beams run different distances, the interferometer must be calibrated in the absence of plasma (numerical simulation number 1). Thus, knowing the phase difference in the absence of plasma, the phase difference due to the plasma is calculated and the electron density as well, from Equation (10). The figure on the right-top represents the calibrated differential phase for all lines of sight, where the phase for N0 = 0 m −3 has been placed equal to zero (calibration). The figure located at the right-bottom shows the true electron density, calculated as the integral on the line of sight of the electron density field, and the electron density calculated by the interferometer, Nine numerical simulations were performed changing the electron density to simulate different plasma conditions, from the absence of the plasma to the flat top condition. The electron density (N 0 ) goes from 0 to 2 × 10 20 m −3 . Figure 4 shows the results of the numerical simulations. The figure on the left represents the ray trajectories from the laser to the detector, where the ray colour shows the calculated phase according to Equation (8). Since the "reference" and the "plasma" beams run different distances, the interferometer must be calibrated in the absence of plasma (numerical simulation number 1). Thus, knowing the phase difference in the absence of plasma, the phase difference due to the plasma is calculated and the electron density as well, from Equation (10). The figure on the right-top represents the calibrated differential phase for all lines of sight, where the phase for N 0 = 0 m −3 has been placed equal to zero (calibration). The figure located at the right-bottom shows the true electron density, calculated as the integral on the line of sight of the electron density field, and the electron density calculated by the interferometer, through the phase difference. The comparison clearly shows that the methodologies implemented to model the interferometer in COMSOL work properly.

Polarimeter Simulations
Analogue numerical simulations were performed to test the polarimeter mod where the plasma current scales accordingly to the electron density. The geometry is s

Polarimeter Simulations
Analogue numerical simulations were performed to test the polarimeter model, where the plasma current scales accordingly to the electron density. The geometry is similar to the interferometer one, where the reference beam is not needed. The three lasers emit beams with a wavelength of 20 µm, and each beam is reflected by a 45 • mirror. Then, they cross a linear polariser ( Figure 5) which sets the polarisation of the probing beam (in the numerical simulation, a linear polarisation with 45 • has been used, since this configuration maximizes the accuracy and the amount of information which can be extracted from the polarimeter, as discussed in the literature [36]). The rays then cross the plasma and their polarisation changes accordingly with Equation (17). The effect of the mirror and linear polariser on the polarisation of the ray are considered implementing the Mueller matrix of each component on the numerical simulation. The left plot of Figure 5 shows the ray trajectories, while the colour of the ray represents the local Cotton-Mouton phase shift. Note that the influences of the mirrors are not represented in the figure, in order to highlight the plasma effect, while for the proper calculation of the plasma parameters from the polarimeter, a calculation of the Mueller matrix (or Jones matrix and others) of the optics must be performed to ensure accurate results. For example, the polarimeter of JET performs a polarisation calibration in the absence of plasma before each plasma discharge. In these numerical simulations, the effects of the optics are calibrated with the first numerical simulation, when there is no plasma. Validation of the Stokes model can be performed by measuring the line integral the electron density with the Cotton-Mouton phase shift, by using the type 1 approxim tion shown in Equation (15). Figure 6 shows the results of this comparison, also demo strating that the polarisation model works properly. Contrariwise, the Faraday rotati also is proportional to the poloidal field (along the propagation direction), making t direct measurement of the electron density difficult. However, being that the Cotto Mouton and the Faraday rotation are linked by the system of Equation (17), a validati on the Cotton-Mouton directly involves the proper propagation of the Faraday rotatio

Thomson Scattering Simulations
The numerical simulations to validate the Thomson scattering model were pe formed changing the geometry according to Figure 7. The light beam, generated by t laser, crosses the plasma. The laser wavelength is 1064 nm, while the peak power is 5× 1 W (note that the TS diagnostic requires a pulsed laser with high peak power). Four diffe The graphs on the right-top and right-bottom show, respectively, the Faraday rotation and the Cotton-Mouton phase shift performed on the nine numerical simulations. The Faraday rotation of the line of sight 1 and 2 are both positive, being the magnetic field along the propagation direction concordant with the light propagation direction, while the line of sight 3 has a negative rotation. Moreover, the Faraday increase is parabolic since it is composed like the product of the electron density and poloidal magnetic field. Contrariwise, the Cotton-Mouton phase shift is positive and increases linearly in each line of sight since it is dominated by the Ω 1 term, which is proportional to the electron density and the toroidal field (the poloidal field influence is negligible in this configuration).
Validation of the Stokes model can be performed by measuring the line integral of the electron density with the Cotton-Mouton phase shift, by using the type 1 approximation shown in Equation (15). Figure 6 shows the results of this comparison, also demonstrating that the polarisation model works properly. Contrariwise, the Faraday rotation also is proportional to the poloidal field (along the propagation direction), making the direct measurement of the electron density difficult. However, being that the Cotton-Mouton and the Faraday rotation are linked by the system of Equation (17) strating that the polarisation model works properly. Contrariwise, the Fa also is proportional to the poloidal field (along the propagation directio direct measurement of the electron density difficult. However, being th Mouton and the Faraday rotation are linked by the system of Equation (17 on the Cotton-Mouton directly involves the proper propagation of the Far

Thomson Scattering Simulations
The numerical simulations to validate the Thomson scattering mo formed changing the geometry according to Figure 7. The light beam, ge laser, crosses the plasma. The laser wavelength is 1064 nm, while the peak W (note that the TS diagnostic requires a pulsed laser with high peak powe ent scattering volumes were placed along the direction of the laser beam a the variables needed for the scattering simulation. For each wavelength, volumes release rays with a specific solid angle aperture. The rays cross th the vacuum vessel, interact with the spherical mirror on the top of the tok rays are collected by the spectrometers, one for each scattering volume.

Thomson Scattering Simulations
The numerical simulations to validate the Thomson scattering model were performed changing the geometry according to Figure 7. The light beam, generated by the laser, crosses the plasma. The laser wavelength is 1064 nm, while the peak power is 5 × 10 8 W (note that the TS diagnostic requires a pulsed laser with high peak power). Four different scattering volumes were placed along the direction of the laser beam and accumulate the variables needed for the scattering simulation. For each wavelength, the scattering volumes release rays with a specific solid angle aperture. The rays cross the plasma, exit the vacuum vessel, interact with the spherical mirror on the top of the tokamak and the rays are collected by the spectrometers, one for each scattering volume. The numerical simulations were conducted changing the N0, T0 and Ip. The density goes from 0 to 20 10 20 m −3 , the electron temperature from 0 to 20 keV and the plasma current from 0 to 15 MA. Figure 8 compares the expected and the simulated spectral density (defined as spectral power divided by the total power, considering the wavelength step equal to 50 nm). For each channel, four different spectra are reported using 0.25, 0.5, 0.75 and 1 times the maximum values of electron density, electron temperature and plasma current. The simulated spectra are perfectly overlapped to the expected spectra at the temperature of each scattering volume, showing that the model works properly. The numerical simulations were conducted changing the N 0 , T 0 and I p . The density goes from 0 to 20 × 10 20 m −3 , the electron temperature from 0 to 20 keV and the plasma current from 0 to 15 MA. Figure 8 compares the expected and the simulated spectral density (defined as spectral power divided by the total power, considering the wavelength step equal to 50 nm). For each channel, four different spectra are reported using 0.25, 0.5, 0.75 and 1 times the maximum values of electron density, electron temperature and plasma current. The simulated spectra are perfectly overlapped to the expected spectra at the temperature of each scattering volume, showing that the model works properly.
Regarding the polarimetric Thomson scattering, the analyses were performed on the same scattering volumes with the same electron, temperature, plasma current and magnetic field. The reduced Stokes vector of the laser beam was placed equal to [1, 0, 0] (linear polarisation with an angle at 0 • ), which is the polarisation normally used in TS diagnostics. In Figure 7, the colour of the ray trajectories shows the first component of the scattered Stokes vector, which decreases as a function of the local temperature. The measurement of the temperature through the polarimetric TS can be performed by using Equation (22).
In the specific case of the numerical simulation, it can be demonstrated that the following equation is valid: Thus, knowing the scattering angle and measuring the first component of the scattered Stokes vector, the temperature can be directly measured. Figure 9 (left) shows the first component of the reduced Stokes vector of the scattered light for the four channels as a function of the electron temperature. Note that the points do not lie on the same line because of the different scattering angle. The other two Stokes components are not shown since they are equal to zero (as expected when the input polarisation has a Stokes vector equal to [1, 0, 0]). Figure 9 (right) shows the expected electron temperature vs. the "measured" temperature by Equation (24)   The numerical simulations were conducted changing the N0, T0 and Ip. The density goes from 0 to 20 10 20 m −3 , the electron temperature from 0 to 20 keV and the plasma current from 0 to 15 MA. Figure 8 compares the expected and the simulated spectral density (defined as spectral power divided by the total power, considering the wavelength step equal to 50 nm). For each channel, four different spectra are reported using 0.25, 0.5, 0.75 and 1 times the maximum values of electron density, electron temperature and plasma current. The simulated spectra are perfectly overlapped to the expected spectra at the temperature of each scattering volume, showing that the model works properly.  because of the different scattering angle. The other two Stokes components are not shown since they are equal to zero (as expected when the input polarisation has a Stokes vector equal to [1, 0, 0]). Figure 9 (right) shows the expected electron temperature vs. the "measured" temperature by Equation (24) using the first component of the reduced Stokes vector. Even in this case, the measured temperature coincides with the expected one, validating the proper functioning of the simulation.

Conclusions
The simulation of diagnostics is a fundamental step in the design process of tokamaks, in order to ensure that the best performances are obtained, and that the degradation due to electromechanical, thermal and nuclear effects is the most limited as possible. These simulations require the use of multiphysics models which can take into account all the effects of degradation and misalignment.
In this work, the modelling of interferometry, polarimetry and Thomson scattering effects have been introduced and implemented on the ray optics module of COMSOL multiphysics. COMSOL was chosen since it allows a very easy interface of the various modules (thermal, mechanics, electromagnetic, etc.), reducing the efforts required to run a multiphysics simulation. However, all the models can be implemented on other platforms.
The validation of the models has been performed by comparing the expected value with the simulated one. All the effects implemented in the model have been tested in an

Conclusions
The simulation of diagnostics is a fundamental step in the design process of tokamaks, in order to ensure that the best performances are obtained, and that the degradation due to electromechanical, thermal and nuclear effects is the most limited as possible. These simulations require the use of multiphysics models which can take into account all the effects of degradation and misalignment.
In this work, the modelling of interferometry, polarimetry and Thomson scattering effects have been introduced and implemented on the ray optics module of COMSOL multiphysics. COMSOL was chosen since it allows a very easy interface of the various modules (thermal, mechanics, electromagnetic, etc.), reducing the efforts required to run a multiphysics simulation. However, all the models can be implemented on other platforms.
The validation of the models has been performed by comparing the expected value with the simulated one. All the effects implemented in the model have been tested in an ITER-like geometry, using a simple circular plasma model (which does not affect the functioning and the validation of the model) with properties comparable with the ones expected in ITER. In the case of the interferometry, three lines of sight have been used. The phase differences between the probing beams and the reference beams, and thus the line integrated electron densities, were calculated. They were compared with the line integral of the electron density of each line of sight, and it was shown that the data are overlapped. In the case of plasma polarimetry, the same probing beams of the interferometer were used, and the simulated Faraday rotation and Cotton-Mouton phase shift were calculated. Then, the line integral density was been calculated using the type 1 approximation of the Cotton-Mouton phase shift. This value was been compared with both the expected and the interferometer quantities, showing perfect accordance and also validating the polarimetric model. The spectral Thomson scattering, since it provides local information, was validated by just comparing simulated spectrum with the expected one, calculated from the local values, while the polarimetric TS was validated calculating the electron temperature from the simulated depolarisation of the scattered radiation and by comparing it with the temperature of the scattering volumes.
This work aimed to analyse the possibility to accurately simulate the optical laserbased diagnostics of thermonuclear plasma and to describe how the model can be implemented in a ray optics code. The validation section clearly shows that the phenomena can be simulated with no errors (between the expected or imposed and the simulated values). In the future, these models may be used to design, test and evaluate the performances of nuclear reactor diagnostics in several conditions, helping in understanding the best setup of the instruments to ensure the best performances and their longest lifetime.