There is a tools within the open-source WhiteboxTools geospatial library called SelectTilesByPolygon that does exactly what you are asking for. It takes an input directory containing many LiDAR files (LAS files) and an input Shapefile polygon and outputs each of the overlapping LiDAR tiles into a separate output directory. The Shapefile can contain multiple polygons.
In the example below, each green square corresponds to a LAS file contained within a directory. The tool extracted each of the 41 LAS files that overlap with the red polygon (a digitized Shapefile) also drawn on the map.
The tool can be run by command-line, using a simple Python GUI (see below), from the QGIS WhiteboxTools plugin, or from a Python script.
Here is the code needed to script the tool (tutorial found here):
from WBT.whitebox_tools import WhiteboxTools
wbt = WhiteboxTools()
las_files_dir = "/path/to/input/files/"
outdir = "/path/to/output/files/"
polygons = "/path/to/input/files/boundary_poly.shp"
wbt.select_tiles_by_polygon(las_files_dir, outdir, polygons)
The only issue that you may have is that WhiteboxTools does not support the compressed LAZ format, just the common LAS file format. The library also contains the ClipLidartoPolygon (for removing points in a LAS file that fall outside a defined region) and ErasePolygonFromLidar (the opposite) tools, which may be of interest for your purpose.
WhiteboxTools can be downloaded from here and a tutorial on processing LiDAR data can be found here. More details are available in the user manual.